Rosen梯度投影法重點講義

1、package XU;import Jama.Matrix;/* * 通用性說明： * 當(dāng)目標(biāo)函數(shù)不同時，程序需修改的地方如下 * 1、函數(shù)getFunction_xy中的f * 2、求梯度的函數(shù)getGradient * 3、線性約束方程組系數(shù)矩陣A * 4、線性約束方程組矩陣b * 5、可行點，初始可行點最好多選擇幾個不同的去求最值，以避免求出的只是區(qū)域極值而不是全域最值 */public class Rosen /實現(xiàn)返回函數(shù)的代數(shù)式，在最優(yōu)化目標(biāo)函數(shù)的表達式變化時，在這個函數(shù)中改即可，避免在進退法和黃金分割法中更改public static double getFunction_xy(

2、double x,double y)double f=0;f=Math.pow(x,2)*y*(4-x-y);return f;/求梯度,注意不同的目標(biāo)函數(shù)，梯度不一樣，此函數(shù)不具有通用性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;/線性約束方程組系數(shù)矩陣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);/以下為各種不需修改的功能函數(shù)/*功能：解決二元非線性函數(shù)在線性約束條件下的極值問題（投影梯度算法） * 注意以下幾點（二元情況下）： * 1、可行點Xi一定為兩個元素的向量 * 2、矩陣A1的大小是變化的，行數(shù)是不定的，但列數(shù)一定為2。A2的情況與A1一樣。b1、b2的行數(shù)是變化的，但其列數(shù)必為1。 * 3、矩陣P的大小一定為2x2 * 4、梯度矩陣的大小一定為2x1 * 5、矩陣d一定為2x1 * 6、矩陣w的行數(shù)與A1相同，列數(shù)為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的行長度（行數(shù)）/int len_A1column_array=A1_array0.length;/A1的列長度（列數(shù)）/int len_A2row_array=A2_array.length;/int len_A2column_array=A2_array0.length;int lenA1_array=dGetlen(A1_array);/lenA1_array0,lenA1_array1分別為行數(shù)和列數(shù)int lenA2_array=dGetlen(A2_array);/lenA2_array0,lenA

7、2_array1分別為行數(shù)和列數(shù)int len_b1_array=b1_array.length;int len_b2_array=b2_array.length;Matrix d=new Matrix(2,1);/創(chuàng)建矩陣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);/創(chuàng)建Matrix類，矩陣中的元素為0.P=Matrix.identity(2,2); /創(chuàng)建單位矩陣P/*System.

8、out.println("矩陣P初始化為：");P.print(2,8); /打印矩陣P*/刪去后i行A1sub_array=getAi(A1_array,i);/*System.out.println("數(shù)組A1sub_array"+i+"為");dPrint(A1sub_array);*/if(k>0)/當(dāng)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(); /矩陣轉(zhuǎn)置./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范數(shù)為所有元素平方和的開方double d_F=d.normF();/System.out.println();/System.out.println("矩陣d"+i+"的F范數(shù)為：");/System.out.println(d_F);/System.out.println();if(d_F<0.00000000000001) /d=0,轉(zhuǎn)第算法中第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(); /矩陣轉(zhuǎn)置./A1sub_T.print(2,8);Matrix w=new Matrix(k,1);/創(chuàng)建矩陣,行數(shù)與A1sub矩陣的行數(shù)相同

13、，列數(shù)為1列w=(A1sub.times(A1sub_T).inverse().times(A1sub).times(Gradient);/System.out.println("矩陣w"+i+"為：");/w.print(lenA1_array0-i,2);/檢測w是否=0,當(dāng)有元素小于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中，最小元素的行標(biāo)值存在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中，最小元素的行標(biāo)值存在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轉(zhuǎn)換為一維向量d_arraydouble d_array=new double2;for(int j=0;j<2;j

18、+)d_arrayj=d.get(j,0);/將矩陣bb轉(zhuǎn)換為一維向量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;/存儲包含極值的的區(qū)間,左右端點分別存放于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目標(biāo)函數(shù)取最小值時的自變量值:xminfc=minHJ(Xi,d_array,0,interval1,0.00000000000001);/Syst

20、em.out.println("目標(biāo)函數(shù)取最小值時的自變量值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;/功能：轉(zhuǎn)置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);/*/打印二位數(shù)組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("數(shù)

24、組為空");System.out.println();/打印一維數(shù)組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();*/矩陣分離函數(shù)(獲得A1)/*函數(shù)說明： * 首先判斷A中的行向量是否滿足式：AX=b，然后將滿足的元素按順序存入一維向量A1_1中， * 然后在將A1_1中的元素按順序存入多維矩陣A1中，最后函數(shù)返回A1，獲得A1矩陣，

25、實現(xiàn)矩陣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為中轉(zhuǎn)存儲矩陣A中滿足式Ai0*Xi0+Ai1*Xi1)=bi的元素的一維矩陣中元素的下標(biāo)值 * m為將中轉(zhuǎn)存儲一維矩陣中的元素存入A1中時的下表計數(shù)變量 * len_h,len_l,分別存儲矩陣的行數(shù)和列數(shù) */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+; /記錄滿足條件的行數(shù)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循環(huán)之后,變量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;/矩陣分離函數(shù)(獲得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+; /記錄滿足條件的行數(shù)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;/矩陣分離函數(shù)(獲得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中的元素的下標(biāo)存在矩陣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;/矩陣分離函數(shù)(獲得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中的元素的下標(biāo)存在矩陣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;/*局限說明：此函數(shù)功能不具有通用性，只針對函數(shù)f=2*x12+2*x22-2*x1*x2-4*x1-6

34、*x2,若要求別的函數(shù)的極值區(qū)間則要 * 在代碼中對f1和f4的表達式進行更改。 * 實現(xiàn)進退法的函數(shù)名為minJT(此處寫的是與Rosen方法配套的函數(shù)，不是單獨的實現(xiàn)進退方法) * 功能:用進退法求解一維函數(shù)的極值區(qū)間 * 參數(shù)說明： * XiRosen方法中迭代的點向量 * dRosen方法中第i次迭代后的d * x0:初始點 * h0:初始步長 * minx:目標(biāo)函數(shù)取包含極值的區(qū)間左端點 * maxx：目標(biāo)函數(shù)取包含極值的區(qū)間右端點 * interval：將左端點和右端點存于數(shù)組interval中 * 關(guān)于進退法詳情參考精通MATLAB最優(yōu)化計算（第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;/* * 實現(xiàn)黃金分割發(fā)函數(shù)名：minHJ * 功能：用黃金分割法求解