Rosen梯度投影法

package XU; import Jama.Matrix;

/* * 通用性说明： * 当目标函数不同时，程序需修改的地方如下 * 1、函数getFunction_xy中的f * 2、求梯度的函数getGradient * 3、线性约束方程组系数矩阵A * 4、线性约束方程组矩阵b * 5、可行点，初始可行点最好多选择几个不同的去求最值，以避免求出的只是区域极值而不是全域最值 */

public class Rosen { //实现返回函数的代数式，在最优化目标函数的表达式变化时，在这个函数中改即可，避免在进退法和黄金分割法中更改 public static double getFunction_xy(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 double[Xi.length]; Gradient[0]=Xi[0]*Xi[1]*(8-3*Xi[0]-2*Xi[1]); Gradient[1]=Math.pow(Xi[0],2)*(4-Xi[0]-2*Xi[1]); /* System.out.println(); System.out.println("梯度为："); for(int i=0;iSystem.out.print(Gradient[i]+","); System.out.println(); */ return Gradient; }

public static void main(String[] args){ double[] Minf=new double[3]; //线性约束方程组系数矩阵A double[][] A={ {-1,-1}, {1,0}, {0,1} }; //线性约束方程组矩阵b double[] b={-6,0,0}; //存储可行点 double[] Xi={2,2}; Minf=getMminf(A,b,Xi); System.out.println(Minf[0]+","+Minf[1]+","+Minf[2]);

} //以下为各种不需修改的功能函数 /*功能：解决二元非线性函数在线性约束条件下的极值问题（投影梯度算法） * 注意以下几点（二元情况下）： * 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[] getMminf(double[][] A,double[] b,double[] Xi){ //存储极值点 double[] X=new double[2];

//精度 double eps=0.00000000000000001;

boolean bConti=true; while(bConti){ //获得矩阵A1,A2,b1,b2 double[][] A1_array=getA1(A,b,Xi); //A1_a为向量p double[][] 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，则为空 //int len_A1row_array=A1_array.length;//A1的行长度（行数） //int len_A1column_array=A1_array[0].length;//A1的列长度（列数） //int len_A2row_array=A2_array.length; //int len_A2column_array=A2_array[0].length; int[] lenA1_array=dGetlen(A1_array);//lenA1_array[0],lenA1_array[1]分别为行数和列数 int[] lenA2_array=dGetlen(A2_array);//lenA2_array[0],lenA2_array[1]分别为行数和列数 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 double[lenA1_array[0]-i][lenA1_array[1]]; int k=lenA1_array[0]-i; Matrix P=new Matrix(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.println("数组A1sub_array"+i+"为"); dPrint(A1sub_array); */ if(k>0){ //当A1不为空

Matrix A1sub=new Matrix(A1sub_array); /* System.out.println("矩阵A1sub"+i+"为"); A1sub.print(2,8); */ Matrix A1sub_T=new Matrix(lenA1_array[1],lenA1_array[0]-i);

A1sub_T=A1sub.transpose(); //矩阵转置. //System.out.println("矩阵A1sub_T"+i+"为："); //A1sub_T.print(lenA1_array[0],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); } 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"+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是否为0 if(k==0){ //A1为空，则找到极值点，算法停止 for(int j=0;j<2;j++) X[j]=Xi[j]; bConti=false; break; } else{

Matrix A1sub=new Matrix(A1sub_array); //System.out.println("矩阵A1sub"+i+"为"); //A1sub.print(2,8); Matrix A1sub_T=new Matrix(lenA1_array[1],lenA1_array[0]-i); A1sub_T=A1sub.transpose(); //矩阵转置. //A1sub_T.print(2,8);

Matrix w=new Matrix(k,1);//创建矩阵,行数与A1sub矩阵的行数相同，列数为1列

w=(((A1sub.times(A1sub_T)).inverse()).times(A1sub)).times(Gradient); //System.out.println("矩阵w"+i+"为："); //w.print(lenA1_array[0]-i,2);