methodOAA.rar

from numpy import * import numpy as np def methodOAA(matrixA,vectorb,vectorc): #步骤1：初始化 alpha=0.999 epsilon=0.001 column=matrixA.shape[1] matrixones=mat(ones((1,column))) bigmatrixones=mat(ones((1,(column+1)))).T #充分大的大M M=900000000000 #带人工变量的增广问题表示成标准形式 vectorc=vstack((vectorc,M)) matrixA1=vectorb-matrixA*matrixones.T matrixA=hstack((matrixA,matrixA1)) #定义单位向量初始可行解,及对应的对角矩阵 solution=np.linspace(1,1,(column+1)) matrixSol=np.mat(np.diag(solution)) solution=np.mat(solution).T I=np.eye(column+1); I=mat(I); #步骤2：计算对偶估计向量 vectorpk=(matrixA*matrixSol*matrixSol*matrixA.T).I*matrixA*matrixSol*matrixSol*vectorc #步骤3：计算即约费用向量 vectorr=vectorc-matrixA.T*vectorpk #步骤4：检查最优性 midvalue=bigmatrixones.T*matrixSol*vectorr print('',midvalue) while vectorr.min()<0 or midvalue>=epsilon: #计算对数壁垒函数目标函数负梯度方向 gradient=-vectorc+0.000000001*1/solution P_k=I-matrixSol*matrixA.T*(matrixA*matrixSol*matrixSol*matrixA.T).I*matrixA*matrixSol #求转移方向 dir=matrixSol*P_k*matrixSol*gradient #判断是否有界 if dir.min()>0: print("原问题无界") break elif dir.min()==0: print('当前解为最优解,solution={}',solution) break else: #计算步长，得到改进解 vectordyk=matrixSol.I*dir minvectordyk=vectordyk[vectordyk<0] distance=(-alpha/minvectordyk).min() solution=solution+distance*dir solution_a=array(solution) solution_a=solution_a.reshape(1,5) #求对角矩阵matrixSol matrixSol=np.diag(solution_a[0]) matrixSol=np.mat(matrixSol) #重复第2步，计算对偶估计向量 vectorpk=(matrixA*matrixSol*matrixSol*matrixA.T).I*matrixA*matrixSol*matrixSol*vectorc #重复第3步，计算即约费用向量 vectorr=vectorc-matrixA.T*vectorpk #重复第四步，检查最优性 midvalue=bigmatrixones.T*matrixSol*vectorr print('目标函数值为{}',vectorc.T*solution) return solution #求解例5.2 A5_2=np.array([[1,-1,1,0],[0,1,0,1]]) A5_2=mat(A5_2) b5_2=mat([[15],[15]]) c5_2=mat([[-2],[1],[0],[0]]) sol=methodOAA(A5_2,b5_2,c5_2); print('最优解为{}',sol)