Rosen梯度投影法

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

2、double x,double y)double f=0;f=Math.pow(x,2)*y*(4-x-y);return f;/求梯度,注意不同的目标函数，梯度不一样，此函数不具有通用性public static double getGradient(double Xi)double Gradient=new doubleXi.length;Gradient0=Xi0*Xi1*(8-3*Xi0-2*Xi1);Gradient1=Math.pow(Xi0,2)*(4-Xi0-2*Xi1);/*System.out.println();System.out.println("梯度为：&

3、quot;);for(int i=0;i<Xi.length;i+)System.out.print(Gradienti+",");System.out.println();*/return Gradient;public static void main(String args)double Minf=new double3;/线性约束方程组系数矩阵Adouble A=-1,-1,1,0,0,1;/线性约束方程组矩阵bdouble b=-6,0,0;/存储可行点double Xi=2,2;Minf=getMminf(A,b,Xi);System.out.printl

4、n(Minf0+","+Minf1+","+Minf2);/以下为各种不需修改的功能函数/*功能：解决二元非线性函数在线性约束条件下的极值问题（投影梯度算法） * 注意以下几点（二元情况下）： * 1、可行点Xi一定为两个元素的向量 * 2、矩阵A1的大小是变化的，行数是不定的，但列数一定为2。A2的情况与A1一样。b1、b2的行数是变化的，但其列数必为1。 * 3、矩阵P的大小一定为2x2 * 4、梯度矩阵的大小一定为2x1 * 5、矩阵d一定为2x1 * 6、矩阵w的行数与A1相同，列数为1列 */public static double getMm

5、inf(double A,double b,double Xi)/存储极值点double X=new double2;/精度double eps=0.00000000000000001;boolean bConti=true;while(bConti)/获得矩阵A1,A2,b1,b2double A1_array=getA1(A,b,Xi); /A1_a为向量pdouble A2_array=getA2(A,b,Xi);double b1_array=getb1(A,b,Xi);double b2_array=getb2(A,b,Xi);/获得矩阵A1,A2,b1,b2的长度，如果长度为0，则

6、为空/int len_A1row_array=A1_array.length;/A1的行长度（行数）/int len_A1column_array=A1_array0.length;/A1的列长度（列数）/int len_A2row_array=A2_array.length;/int len_A2column_array=A2_array0.length;int lenA1_array=dGetlen(A1_array);/lenA1_array0,lenA1_array1分别为行数和列数int lenA2_array=dGetlen(A2_array);/lenA2_array0,lenA

7、2_array1分别为行数和列数int len_b1_array=b1_array.length;int len_b2_array=b2_array.length;Matrix d=new Matrix(2,1);/创建矩阵d,2行1列int i=0;while(true)double A1sub_array=new doublelenA1_array0-ilenA1_array1;int k=lenA1_array0-i;Matrix P=new Matrix(2,2);/创建Matrix类，矩阵中的元素为0.P=Matrix.identity(2,2); /创建单位矩阵P/*System.

8、out.println("矩阵P初始化为：");P.print(2,8); /打印矩阵P*/删去后i行A1sub_array=getAi(A1_array,i);/*System.out.println("数组A1sub_array"+i+"为");dPrint(A1sub_array);*/if(k>0)/当A1不为空Matrix A1sub=new Matrix(A1sub_array);/*System.out.println("矩阵A1sub"+i+"为");A1sub.print

9、(2,8);*/Matrix A1sub_T=new Matrix(lenA1_array1,lenA1_array0-i);A1sub_T=A1sub.transpose(); /矩阵转置./System.out.println("矩阵A1sub_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("矩阵

10、P"+i+"为：");/P.print(2,8);double Gradient_arraay=getGradient(Xi);Matrix Gradient=new Matrix(Gradient_arraay,2);/获得2x1梯度矩阵/System.out.println();/System.out.println("梯度矩阵"+i+"为：");/Gradient.print(1,2);d=(P.times(-1).times(Gradient);/计算dk/System.out.println("矩阵d&q

11、uot;+i+"为：");/d.print(1,8);/检测d是否为0,F范数为所有元素平方和的开方double d_F=d.normF();/System.out.println();/System.out.println("矩阵d"+i+"的F范数为：");/System.out.println(d_F);/System.out.println();if(d_F<0.00000000000001) /d=0,转第算法中第5步 ，检测d是否为0if(k=0)/A1为空，则找到极值点，算法停止for(int j=0;j<2

12、;j+)Xj=Xij;bConti=false;break;elseMatrix A1sub=new Matrix(A1sub_array);/System.out.println("矩阵A1sub"+i+"为");/A1sub.print(2,8);Matrix A1sub_T=new Matrix(lenA1_array1,lenA1_array0-i);A1sub_T=A1sub.transpose(); /矩阵转置./A1sub_T.print(2,8);Matrix w=new Matrix(k,1);/创建矩阵,行数与A1sub矩阵的行数相同

13、，列数为1列w=(A1sub.times(A1sub_T).inverse().times(A1sub).times(Gradient);/System.out.println("矩阵w"+i+"为：");/w.print(lenA1_array0-i,2);/检测w是否=0,当有元素小于0时，w_0为falseboolean w_0=true;for(int j=0;j<k;j+)if(w.get(j,0)<0)w_0=false;/System.out.println("检测w"+i+"是否=0的布尔变量值（

14、true表明=0，false表明有元素<0）为：");/System.out.println(w_0);if(w_0) /w>=0,则找到极值点，算法停止for(int j=0;j<2;j+)Xi=Xii;bConti=false;break;else/找出w中最小的元素，并存放在w_array0中，最小元素的行标值存在index中double w_array=new doublek;int index=-1;for(int j=0;j<k;j+)w_arrayj=w.get(j,0);for(int j=0;j<k;j+)if(w_arrayj<

15、=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/System.out.println("哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈");break;/以下部分处理算法中的一维部分Matrix bb=new Matrix(len_b2_arr

16、ay,1);Matrix A2=new Matrix(A2_array);Matrix b2=new Matrix(b2_array,len_b2_array);Matrix Xi_M=new Matrix(Xi,2);Matrix dd=new Matrix(lenA2_array0,1);bb=b2.minus(A2.times(Xi_M);/计算b'=b2-A2*Xdd=A2.times(d);/计算d'=A2*d/bb.print(2,2);/dd.print(1,2);/找出dd中最小的元素，并存放在dd_array0中，最小元素的行标值存在index1中double

17、 dd_array_min=0;double dd_array=new doublelenA2_array0;int index1=-1;for(int j=0;j<lenA2_array0;j+)dd_arrayj=dd.get(j,0);for(int j=0;j<lenA2_array0;j+)if(dd_arrayj<=dd_array0)dd_array_min=dd_arrayj;index1=j;/sPrint(dd_array);/将矩阵d转换为一维向量d_arraydouble d_array=new double2;for(int j=0;j<2;j

18、+)d_arrayj=d.get(j,0);/将矩阵bb转换为一维向量bb_arraydouble bb_array=new doublelen_b2_array;for(int j=0;j<len_b2_array;j+)bb_arrayj=bb.get(j,0);/sPrint(bb_array);double interval=new double2;/存储包含极值的的区间,左右端点分别存放于interval0、interval1中if(dd_array_min>=0)interval=minJT(Xi,d_array,0,0.000001);elseinterval1=10

19、0000;for(int j=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);double minfc=new double2;/minfc0目标函数取最小值时的自变量值:xminfc=minHJ(Xi,d_array,0,interval1,0.00000000000001);/Syst

20、em.out.println("目标函数取最小值时的自变量值x:");/System.out.println(minfc0);double tol=d.times(minfc0).normF();if(tol<eps)X=Xi;break;for(int j=0;j<2;j+)Xij=Xij+minfc0*d_arrayj;/System.out.println("Xi为");/sPrint(Xi);/double minf=2*Math.pow(X0,2)+2*Math.pow(X1,2)-2*X0*X1-4*X0-6*X1;double

21、minf=getFunction_xy(X0,X1);/*System.out.println();for(int j=0;j<2;j+)System.out.print(Xj+",");System.out.println();System.out.println("minf="+minf);*/double Minf=new double3;for(int j=0;j<2;j+)Minfj=Xj;Minf2=minf;return Minf;/功能：转置public static double getA_T(double A1)int m

22、=A1.length;int n=A10.length;double A1_T=new doublenm;for (int i = 0; i<n;i+) for (int j=0; j<m;j+) A1_Tij = A1ji; return A1_T;/功能：去掉A的后h行public static double getAi(double A,int h)int m=A.length;if(m!=0)int n=A0.length;double A_sub=new doublem-hn;for(int i=0;i<m-h;i+)for(int j=0;j<n;j+)A_

23、subij=Aij;return A_sub;elsedouble Kong=new double00;return(Kong);/*/打印二位数组public static void dPrint(double N)int m=N.length;if(m!=0)int n=N0.length;for(int i=0;i<m;i+)for(int j=0;j<n;j+)System.out.print(Nij+" ");System.out.println();System.out.println();elseSystem.out.println("数

24、组为空");System.out.println();/打印一维数组public static void sPrint(double N)int l=N.length;for(int i=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矩阵，

25、实现矩阵A分解。 */public static double getA1(double A,double b,double Xi)int hang=0,num=0; int len_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*le

26、n_l;/*System.out.println("len_h="+len_h);System.out.println("len_l="+len_l);System.out.println("len="+len);*/double A1_1=new doublelen;/声明A1_1，其大小和A一样/将A中满足条件的元素按顺序存入一维向量A1_1中for(int i=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi<0.000000000000000000001)hang+; /记录满足条件的行数fo

27、r(int j=0;j<len_l;j+)A1_1num+=Aij;/System.out.println("num="+num);/System.out.println();/打印A1_1/System.out.println("向量A1_1为：");/sPrint(A1_1);/将A1_1中的元素按顺序放入A1中double A1=new doublehanglen_l;/注意，这句开辟A1空间的语句，一定要放在判断语句块for循环之后,变量hang才会有意义。num=0;for(int i=0;i<hang;i+)for(int j=0

28、;j<2;j+)A1ij=A1_1num+;/打印矩阵A1/System.out.println("向量A1为：");/dPrint(A1);return A1;/矩阵分离函数(获得A2)public static double getA2(double A,double b,double Xi)int hang=0,num=0; int len_h,len_l,len;len_h=A.length;len_l=A0.length;len=len_h*len_l;double A2_1=new doublelen;/将A中不满足条件的元素按顺序存入一维向量A2_1中f

29、or(int i=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi>=0.000000000000000000001)hang+; /记录满足条件的行数for(int j=0;j<len_l;j+)A2_1num+=Aij;/打印A2_1/System.out.println("向量A2_1为：");/sPrint(A2_1);/将A2_1中的元素按顺序放入A2中double A2=new doublehanglen_l;num=0;for(int i=0;i<hang;i+)for(int j=0;j<2;j+)A2ij

30、=A2_1num+;/打印矩阵A2/System.out.println("向量A2为：");/dPrint(A2);return A2;/矩阵分离函数(获得b1)public static double getb1(double A,double b,double Xi)int len_h=A.length;int len_b=b.length;int index=new intlen_b;int num=0;/将满足条件的矩阵b中的元素的下标存在矩阵index中for(int i=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi<0.0

31、00000000000000000001)indexnum+=i;/将b中满足条件的元素存入矩阵b1中double b1=new doublenum;for(int i=0;i<num;i+)b1i=bindexi;/打印b1/System.out.println("向量b1为：");/sPrint(b1);return b1;/矩阵分离函数(获得b2)public static double getb2(double A,double b,double Xi)int len_h=A.length;int len_b=b.length;int index=new in

32、tlen_b;int num=0;/将满足条件的矩阵b中的元素的下标存在矩阵index中for(int i=0;i<len_h;i+)if(Ai0*Xi0+Ai1*Xi1)-bi>=0.000000000000000000001)indexnum+=i;/将b中满足条件的元素存入矩阵b1中double b2=new doublenum;for(int i=0;i<num;i+)b2i=bindexi;/打印b1/System.out.println("向量b2为：");/sPrint(b2);return b2;/矩阵i、j行对换public static

33、 double rowChange(double N,int i,int j)int m=N.length;int n=N0.length;double N_change=new doublemn;for(int l=0;l<m;l+)for(int r=0;r<n;r+)N_changelr=Nlr;for(int r=0;r<n;r+)N_changejr=Nir;for(int r=0;r<n;r+)N_changeir=Njr;return N_change;/*局限说明：此函数功能不具有通用性，只针对函数f=2*x12+2*x22-2*x1*x2-4*x1-6

34、*x2,若要求别的函数的极值区间则要 * 在代码中对f1和f4的表达式进行更改。 * 实现进退法的函数名为minJT(此处写的是与Rosen方法配套的函数，不是单独的实现进退方法) * 功能:用进退法求解一维函数的极值区间 * 参数说明： * XiRosen方法中迭代的点向量 * dRosen方法中第i次迭代后的d * x0:初始点 * h0:初始步长 * minx:目标函数取包含极值的区间左端点 * maxx：目标函数取包含极值的区间右端点 * interval：将左端点和右端点存于数组interval中 * 关于进退法详情参考精通MATLAB最优化计算（第2版）.龚纯等 */public

35、static double minJT(double Xi,double d_array,double x0,double h0)double x1=0,x2=0,x3=0,x4=0,h=0,f1=0,f2=0,f3=0,f4=0,minx=0,maxx=0;double interval=new double2;x1=x0;int k=0;h=h0;while(true)x4=x1+h;k=k+1;/求f4double X_4=new double2;X_40=Xi0+d_array0*x4;X_41=Xi1+d_array1*x4;/f4=2*X_40*X_40+2*X_41*X_41-2

36、*X_40*X_41-4*X_40-6*X_41;f4=getFunction_xy(X_40,X_41);/求f1double X_1=new double2;X_10=Xi0+d_array0*x1;X_11=Xi1+d_array1*x1;/f1=2*X_10*X_10+2*X_11*X_11-2*X_10*X_11-4*X_10-6*X_11;f1=getFunction_xy(X_10,X_11);if(f4<f1)x2=x1;x1=x4;f2=f1;f1=f4;h=2*h;elseif(k=1)h=-h;x2=x4;f2=f4;elsex3=x2;x2=x1;x1=x4;break;if(x1>x3)minx=x3;elseminx=x1;maxx=x1+x3-minx;interval0=minx;interval1=maxx;return interval;/* * 实现黄金分割发函数名：minHJ * 功能：用黄金分割法求解