Zoutendijk算法python实现

Zoutendijk 可行方向法通过下述问题求搜索方向：

算法流程为：

#目标函数def targetfun(x):    f_d = fgrend(x)    return numpy.dot(f_d,x)#梯度函数def fgrend(x):    f_d = numpy.array([2*x[0]+x[1]-2*x[2]-4,4*x[1]+x[0]+x[2]-6,6*x[2]-2*x[0]+x[1]])    return f_d#Zoutendijk算法def zoutendik(k,x,A,b):    A_effect = zeros(3)    b_effect = zeros(1)    A_effect2 = zeros(3)    b_effect2 = zeros(1)    for i in range(len(A)):        if numpy.dot(A[i],x)==b[i]:            if len(b_effect)==1:                A_effect = A[i]                b_effect = b[i]            else:                A_effect = numpy.row_stack((A_effect,A[i]))                b_effect = numpy.row_stack((b_effect,b[i]))        else:            if len(b_effect2)==1:                A_effect2 = A[i]                b_effect2 = b[i]            else:                A_effect2 = numpy.row_stack((A_effect2,A[i]))                b_effect2 = numpy.row_stack((b_effect2,b[i]))    # 先求搜索方向d    f_d = fgrend(x)    Q =  lambda d:  numpy.dot(f_d , d)    cons = ({'type': 'ineq', 'fun': lambda x:  numpy.dot(A_effect[0], x)},            {'type': 'ineq', 'fun': lambda x:  numpy.dot(A_effect[1], x)},            {'type': 'ineq', 'fun': lambda x:  numpy.dot(A_effect[2], x)},            )    bnds = ((-1, 1), (-1, 1), (-1, 1))    res = minimize(Q, (0, 0,0),   bounds=bnds, constraints=cons)    dvalue = res['x']    z = numpy.dot(f_d,dvalue)    #再求步长lambda     lamb =  lambda l:  (x+ dvalue*l)[0]**2+2*(x+dvalue *l)[1]**2+3*(x+dvalue *l)[2]**2+(x+dvalue *l)[0]* (x+dvalue*l)[1]-2*(x+dvalue *l)[0]*(x+dvalue *l)[2]+(x+dvalue *l)[1]*(x+dvalue *l)[2]-4*(x+dvalue *l)[0]-6*(x+dvalue *l)[1]    des = numpy.dot(A_effect2,dvalue)    bes = b_effect2-numpy.dot(A_effect2,x)    lambdaset = []    des = numpy.array(des)    bes = numpy.array(bes)    mlambda = bes/des    #for i in  range(1):    if des<0:        #lambdaset.append(float(bes[i])/float(des[i]))#一般情形下需要计算比值        max_lambda = mlambda    else:        max_lambda = 1000000    cons = ({'type': 'ineq', 'fun': lambda x: x}             )    bnds = ((0, max_lambda),)#注意，此处如果不加','长度为2    xinit = numpy.array([0])    reslam = minimize(lamb, xinit,   bounds=bnds, constraints=cons)    lambdavalue = reslam['x']    f  = targetfun(x)    print "第%s轮迭代：x=%s,d=%s,lambda=%s，z=%s,f=%s"%(k,x,dvalue,lambdavalue,z,f)    return (z,lambdavalue,dvalue)def controller():    j = 0    x0 = numpy.array([0,0,0])    f_d = fgrend(x0)    A = numpy.array([[-1,-2,-1],[1,0,0],[0,1,0],[0,0,1]])    b = numpy.array([[-4],[0],[0],[0]])    x = x0    zou = zoutendik(j,x0,A,b)    z,lambd,d = zou    for i in range(1,500):        x = x+ lambd*d        zu = zoutendik(i,x,A,b)        z = zu[0]        if z:            zu1 = zoutendik(i,x,A,b)            z,lambd,d = zu1