Rosen梯度投影法

1、packageXU;importJama.Matrix;/*通用性说明：*当目标函数不同时，程序需修改的地方如下*1、函数getFunction_xy中的f*2、求梯度的函数getGradient*3、线性约束方程组系数矩阵A*4、线性约束方程组矩阵b*5、可行点，初始可行点最好多选择几个不同的去求最值，以避免求出的只是区域极值而不是全域最值*/publicclassRosen/实现返回函数的代数式，在最优化目标函数的表达式变化时，在这个函数中改即可，避免在进退法和黄金分割法中更改publicstaticdoublegetFunction_xy(doublex,doubley)doublef=

2、0;f=Math.pow(x,2)*y*(4-x-y);returnf;/求梯度,注意不同的目标函数，梯度不一样，此函数不具有通用性publicstaticdoublegetGradient(doubleXi)doubleGradient=newdoubleXi.length;Gradient0=Xi0*Xi1*(8-3*Xi0-2*Xi1);Gradient1=Math.pow(Xi0,2)*(4-Xi0-2*Xi1);/*System.out.println();System.out.printlnC梯度为：”)；for(inti=0;i<Xi.length;i+)System.ou

3、t.print(Gradienti+",")；System.out.println()；*/returnGradient；publicstaticvoidmain(Stringargs)doubleMinf=newdouble3；线性约束方程组系数矩阵AdoubleA=-1,-1,1,0,0,1;线性约束方程组矩阵bdoubleb=-6,0,0;/存储可行点doubleXi=2,2;Minf=getMminf(A,b,Xi);System.out.println(Minf0+","+Minf1+","+Minf2);/以下为各种不需修

4、改的功能函数/*功能：解决二元非线性函数在线性约束条件下的极值问题(投影梯度算法)*注意以下几点(二元情况下)：*1、可行点Xi一定为两个元素的向量* 2、矩阵A1的大小是变化的，行数是不定的，但列数一定为2。A2的情况与A1一样。bl、b2的行数是变化的，但其列数必为1。* 3、矩阵P的大小一定为2x2* 4、梯度矩阵的大小一定为2x1* 5、矩阵d一定为2x1* 6、矩阵w的行数与A1相同，列数为1列*/publicstaticdoublegetMminf(doubleA,doubleb,doubleXi)/存储极值点doubleX=newdouble2;/精度doubleeps=0.00

5、000000000000001;booleanbConti=true;while(bConti)/获得矩阵A1,A2,b1,b2doubleA1_array=getA1(A,b,Xi);A1_a为向量pdoubleA2_array=getA2(A,b,Xi);doubleb1_array=getb1(A,b,Xi);doubleb2_array=getb2(A,b,Xi);获得矩阵A1,A2,b1,b2的长度，如果长度为0则为空intlen_A1row_array=A1_array.length;/A1的行长度(行数)intlen_A1column_array=A1_array0.length

6、;/A1的列长度(列数)/intlen_A2row_array=A2_array.length;/intlen_A2column_array=A2_array0.length;intlenA1_array=dGetlen(A1_array);/lenA1_array0,lenA1_array1分别为行数和列数intlenA2_array=dGetlen(A2_array);/lenA2_array0,lenA2_array1分别为行数和列数intlen_b1_array=b1_array.length;intlen_b2_array=b2_array.length;Matrixd=newMat

7、rix(2,1);仓U建矩阵d,2行1列inti=0;while(true)doubleA1sub_array=newdoublelenA1_array0-ilenA1_array1;intk=lenA1_array0-i;MatrixP=newMatrix(2,2);仓建Matrix类，矩阵中的元素为0.P=Matrix.identity(2,2);/仓建单位矩阵P/*System.out.println('矩阵P初始化为：")；P.print(2,8);/打印矩阵P*/删去后i行A1sub_array=getAi(A1_array,i);/*System.out.prin

8、tlnC数组A1sub_array"+i+"为")；dPrint(A1sub_array);*/if(k>0)当A1不为空MatrixA1sub=newMatrix(A1sub_array);/*System.out.printlnC矩阵A1sub"+i+"为")；A1sub.print(2,8);*/MatrixA1sub_T=newMatrix(lenA1_array1,lenA1_array0-i);A1sub_T=A1sub.transpose();/矩阵转置.System.out.println('矩阵A1su

9、b_T"+i+"为：”)；/A1sub_T.print(lenA1_array0,2);P.minusEquals(A1sub_T.times(A1sub.times(A1sub_T).inverse().times(A1sub);/计十算P=I-(A1).System.out.println('矩阵P"+i+"为：”)；/P.print(2,8);doubleGradient_arraay=getGradient(Xi);MatrixGradient=newMatrix(Gradient_arraay,2);获得2x1梯度矩阵/System.o

10、ut.println();/System.out.printlnC梯度矩阵"+i+"为：”)；/Gradient.print(1,2);d=(P.times(-1).times(Gradient);计算dkSystem.out.println('矩阵d"+i+"为：”)；/d.print(1,8);检测d是否为0,F范数为所有元素平方和的开方doubled_F=d.normF();/System.out.println();System.out.println('矩阵d"+i+"的F范数为：”)；/System.out

11、.println(d_F);/System.out.println();if(d_F<O.OOOOOOOOOOOOO1)d=0，转第算法中第5步，检测d是否为0if(k=O)/A1为空，则找到极值点，算法停止for(intj=O;j<2;j+)Xj=Xij;bConti=false;break;elseMatrixA1sub=newMatrix(A1sub_array);System.out.println('矩阵A1sub"+i+"为”)；/A1sub.print(2,8);MatrixA1sub_T=newMatrix(lenA1_array1,le

12、nA1_arrayO-i)；A1sub_T=A1sub.transpose()；/矩阵转置./A1sub_T.print(2,8)；Matrixw=newMatrix(k,1);创建矩阵，行数与A1sub矩阵的行数相同，列数为1列w=(A1sub.times(A1sub_T).inverse().times(A1sub).times(Gradient)；System.out.println('矩阵w"+i+"为：”)；/w.print(lenA1_arrayO-i,2)；检测w是否=0,当有元素小于0时，w_0为falsebooleanw_0=true;for(in

13、tj=0;j<k;j+)if(w.get(j,0)<0)w_0=false;System.out.println('检测w"+i+"是否=0的布尔变量值(true表明=0,false表明有元素<0)为：");/System.out.println(w_0);if(w_0)w>=0,则找到极值点，算法停止for(intj=0;j<2;j+)Xi=Xii;bConti=false;break;else找出w中最小的元素，并存放在w_array0冲,最小元素的行标值存在index中doublew_array=newdoublek;i

14、ntindex=-1;for(intj=0;j<k;j+)w_arrayj=w.get(j,0);for(intj=0;j<k;j+)if(w_arrayj<=w_array0)w_array0=w_arrayj;index=j;将A1_array的j行与最后一行对换A1_array=rowChange(A1sub_array,index,k-1);/System.out.println();System.out.println("A1_array的j行与最后一行对换后为");/dPrint(A1_array);i+;else值存在index1中别存放于i

15、nterval0、/System.out.printlnC'哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈")；break;/以下部分处理算法中的一维部分Matrixbb=newMatrix(len_b2_array,1)；MatrixA2=newMatrix(A2_array)；Matrixb2=newMatrix(b2_array,len_b2_array)；MatrixXi_M=newMatrix(Xi,2)；Matrixdd=newMatrix(lenA2_array0,1)；bb=b2.minus(A2.times(Xi_M)；/计算b'=b2-A2*Xdd=A2.time

16、s(d)；/计算d'=A2*d/bb.print(2,2)；/dd.print(1,2)；找出dd中最小的元素，并存放在dd_array0冲，最小元素的行标doubledd_array_min=0；doubledd_array=newdoublelenA2_array0；intindex1=-1；for(intj=0；j<lenA2_array0；j+)dd_arrayj=dd.get(j,0)；for(intj=0；j<lenA2_array0；j+)if(dd_arrayj<=dd_array0)dd_array_min=dd_arrayj；index1=j；/s

17、Print(dd_array)；/将矩阵d转换为一维向量d_arraydoubled_array=newdouble2；for(intj=0；j<2；j+)d_arrayj=d.get(j,0)；/将矩阵bb转换为一维向量bb_arraydoublebb_array=newdoublelen_b2_array；for(intj=0；j<len_b2_array；j+)bb_arrayj=bb.get(j,0)；/sPrint(bb_array)；doubleinterval=newdouble2;存储包含极值的的区间，左右端点分interval中if(dd_array_min>

18、;=0)interval=minJT(Xi,d_array,0,0.000001);elseinterval1=100000;for(intj=0;j<dd_array.length;j+)if(dd_arrayj<0)if(bb_arrayj/dd_arrayj<interval1)interval1=bb_arrayj/dd_arrayj;System.out.println('右端点为")；/System.out.println(interval1);doubleminfc=newdouble2;/minfc0目标函数取最小值时的自变量值:xminfc

19、=minHJ(Xi,d_array,0,interval1,0.00000000000001)；/System.out.printlnC'目标函数取最小值时的自变量值x:");/System.out.println(minfc0);doubletol=d.times(minfc0).normF();if(tol<eps)X=Xi;break;for(intj=0;j<2;j+)Xij=Xij+minfc0*d_arrayj;/System.out.println("Xi为");/sPrint(Xi);/doubleminf=2*Math.pow

20、(X0,2)+2*Math.pow(X1,2)-2*X0*X1-4*X0-6*X1;doubleminf=getFunction_xy(X0,X1);/*System.out.println();for(intj=0;j<2;j+)System.out.print(Xj+",");System.out.println();System.out.println("minf="+minf);*/doubleMinf=newdouble3;for(intj=0;j<2;j+)Minfj=Xj;Minf2=minf;returnMinf;/功能：转置

21、publicstaticdoublegetA_T(doubleA1)intm=A1.length;intn=A10.length;doubleA1_T=newdoublenm;for(inti=0;i<n;i+)for(intj=0;j<m;j+)A1_Tij=A1ji;returnA1_T;/功能：去掉A的后h行publicstaticdoublegetAi(doubleA,inth)intm=A.length;if(m!=0)intn=A0.length;doubleA_sub=newdoublem-hn;for(inti=0;i<m-h;i+)for(intj=0;j&

22、lt;n;j+)A_subij=Aij;returnA_sub;elsedoubleKong=newdouble00;return(Kong);/*/打印二位数组publicstaticvoiddPrint(doubleN)intm=N.length;if(m!=0)intn=N0.length;for(inti=0;i<m;i+)for(intj=0;j<n;j+)System.out.print(Nij+"");System.out.println();System.out.println();elseSystem.out.printlnC数组为空"

23、;)；System.out.println();/打印一维数组publicstaticvoidsPrint(doubleN)intl=N.length；for(inti=0；i<l；i+)System.out.print(Ni+"")；System.out.println()；System.out.println()；*/矩阵分离函数(获得A1)/*函数说明：*首先判断A中的行向量是否满足式：AX=b，然后将满足的元素按顺序存入一维向量A1_1中，*然后在将A1_1中的元素按顺序存入多维矩阵A1中，最后函数返回A1，获得A1矩阵，实现矩阵A分解。*/publicsta

24、ticdoublegetA1(doubleA,doubleb,doubleXi)inthang=0,num=0；intlen_h,len_l,len；/*变量hang为确定A1有几行，hang>0，说明A1不为空*变量num为中转存储矩阵A中满足式Ai0*Xi0+Ai1*Xi1)=bi的元素的一维矩阵中元素的下标值* m为将中转存储一维矩阵中的元素存入A1中时的下表计数变量* len_h,len_l,分别存储矩阵的行数和列数*/len_h=A.length;len_l=A0.length;len=len_h*len_l;/*System.out.println("len_h=&

25、quot;+len_h);System.out.println("len_l="+len_l);System.out.println("len="+len);*/doubleA1_1=newdoublelen;声明A1_1,其大小和A一样将A中满足条件的元素按顺序存入一维向量A1_1中for(inti=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi<0.000000000000000000001)hang+;/记录满足条件的行数for(intj=0;j<len_l;j+)A1_1num+=Aij;/System.

26、out.println("num="+num);/System.out.println();/打印A1_1System.out.println("向量A1_1为：")；/sPrint(A1_1);将A1_1中的元素按顺序放入A1中doubleA1=newdoublehanglen_l;/注意，这句开辟A1空间的语句，一定要放在判断语句块for循环之后，变量hang才会有意义。num=0;for(inti=0;i<hang;i+)for(intj=0;j<2;j+)A1ij=A1_1num+;/打印矩阵A1System.out.println(

27、"向量A1为：")；/dPrint(A1);returnA1;/矩阵分离函数(获得A2)publicstaticdoublegetA2(doubleA,doubleb,doubleXi)inthang=0,num=0;intlen_h,len_l,len;len_h=A.length;len_l=A0.length;len=len_h*len_l;doubleA2_1=newdoublelen;将A中不满足条件的元素按顺序存入一维向量A2_1中for(inti=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi>=0.000000000000

28、000000001)hang+;/记录满足条件的行数for(intj=0;j<len_l;j+)A2_1num+=Aij;/打印A2_1System.out.println("向量A2_1为：")；/sPrint(A2_1);将A2_1中的元素按顺序放入A2中doubleA2=newdoublehanglen_l；num=0；for(inti=0；i<hang；i+)for(intj=0；j<2；j+)A2ij=A2_1num+；/打印矩阵A2System.out.println("向量A2为：")；/dPrint(A2)；return

29、A2；/矩阵分离函数(获得b1)publicstaticdoublegetb1(doubleA,doubleb,doubleXi)intlen_h=A.length；intlen_b=b.length；intindex=newintlen_b；intnum=0；将满足条件的矩阵b中的元素的下标存在矩阵index中for(inti=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi<0.000000000000000000001)indexnum+=i;将b中满足条件的元素存入矩阵bl中doubleb1=newdoublenum;for(inti=0;i<n

30、um;i+)b1i=bindexi;打印blSystem.out.println("向量bl为：")；/sPrint(b1);returnb1；/矩阵分离函数(获得b2)publicstaticdoublegetb2(doubleA,doubleb,doubleXi)intlen_h=A.length；intlen_b=b.length；intindex=newintlen_b；intnum=0；将满足条件的矩阵b中的元素的下标存在矩阵index中for(inti=0；i<len_h；i+)if(Ai0*Xi0+Ai1*Xi1)-bi>=0.0000000000

31、00000000001)indexnum+=i；将b中满足条件的元素存入矩阵bl中doubleb2=newdoublenum；for(inti=0；i<num；i+)b2i=bindexi；/打印b1System.out.println("向量b2为：")；/sPrint(b2)；returnb2；矩阵i、j行对换publicstaticdoublerowChange(doubleN,inti,intj)intm=N.length;intn=N0.length;doubleN_change=newdoublemn;for(intl=0;l<m;l+)for(in

32、tr=0;r<n;r+)N_changelr=Nlr;for(intr=0;r<n;r+)N_changejr=Nir;for(intr=0;r<n;r+)N_changeir=Njr;returnN_change;/*局限说明：此函数功能不具有通用性，只针对函数f=2*x"2+2*x2A2-2*x1*x2-4*x1-6*x2，若要求别的函数的极值区间则要* 在代码中对f1和f4的表达式进行更改。*实现进退法的函数名为minJT（此处写的是与Rosen方法配套的函数，不是单独的实现进退方法）* 功能:用进退法求解一维函数的极值区间*参数说明：*XiRosen方法中迭代的点向量*dRosen方法中第i次迭代后的d*x0:初始点*h0:初始步长*minx:目标函数取包含极值的区间左端点*maxx：目标函数取包含极值的区间右端点*interval：将左端点和右端点存于数组interval中*关于进退法详情参考精通MATLAB最优化计算(第2版)龚纯等*/publicstaticdoubleminJT(doubleXi,doubled_array,doublex0,double