最优化方法课程实验报告

项目一一维搜索算法(一)［实验目的］编写加步探索法、对分法、Newton法的程序。［实验准备］.掌握一维收搜索中搜索区间的加步探索法的思想及迭代步骤;2.掌握对分法的思想及迭代步骤；3.掌握Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：.用加步探索法确定一维最优化问题min9(t)=13—2t+1t>0的搜索区间，要求选取10=°’h0=1,以=2.加步探索法算法的计算步骤：. ..A占tG［0,+8)(或tGO［0，t］)斗曾甲=9(t)外山河仙牛Kh〉0⑴选取初女口点0 0 max，计算0 0.给出初女口步长h0〉0,加步系数以I，令k=0。比较目标函数值.令L=tk+hk，计算9k+i=9(L)，若9k+i<9k，转(3)，否则转(4)。 + +k=k+1加大探索步长.令hk+i以hk，同时，令 k,\tk+i, ，转(2)。⑷反向探索.若k=0，转换探索方向，令hk=—K’t=tk+i，转(2)。否则，停止迭代，aa=min{t,t},b=max{t,t}加步探索法算法的计算框图程序清单加步探索法算法程序见附录1实验结果运行结果为：■"D:\ProgrsmFiles ero^o+tVisualStudio\MyPr■:该问题的根的搜索空间是'[0^1Pressanykeytocontinue—2.用对分法求解min中(t)=t(t+2)，已知初始单谷区间"Ml=［—3,5］，要求按精度'=0.3，£=0.001分别计算.对分法迭代的计算步骤：确定初始搜索区间［a,b］,要求中'(a)<0中'(b)>0。,计算［a,b］的中点c=-(a+b).若甲'(c)<0，则a=c，转(4);若甲'(c)=0，则t*=c，转(5)；若甲'(c)>0，则b=c转⑷.1⑷若Ia—b1<£，则t*=5(a+b)，转(5)；否则转(2).(5)打印，*，结束对分法的计算框图程序清单对分法程序见附录2实验结果运行结果为：LI'D:\PnogramFiles(x86AMicr-csoftViiuaIStudic\MyProject3解的结果是:TPressanykeytocontinue3.用Newton法求解min中(t)=13—2t+1,已知初始单谷区间"Ml=[0，1],要求精度£=°-01.Newton法的计算步骤⑴确定初始搜索区间[a,b]，要求甲'(a)<0,甲'(b)>0⑵选定10(3)计算t='。一甲'("甲"(t°)⑷若It—101>8，则t0=t，转⑶；否则转(5).(5)打印t,(5)打印t,中(t)结束.Newton法的计算框图程序清单Newton法程序见附录3实验结果运行结果为：项目二一维搜索算法（二）［实验目的］编写黄金分割法、抛物线插值法的程序。［实验准备］1.掌握黄金分割法的思想及迭代步骤；2.掌握抛物线插值法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用黄金分割法求解

min中(t)=t(t+2),已知初始单谷区间"Ml=[—3,5]，要求精度£=。.°°1.黄金分割法迭代步骤：确定中(t)的初始搜索区间[a,b].计算12=a+0.382(b-a)计算t1=a+0.618(b-a)若11-1|<e，则打印t*=，结束；否则转(5).TOC\o"1-5"\h\z1 2 2判别是否满足VV修：若满足，则置a=t,t=t，^*4 4 44然后转(3)；否则，置b=t,t=t,V=V，t=以+。(b—a),V=V(t)11 2 1 22 2 2然后转(4).黄金分割法的计算框图：程序清单黄金分割法程序见附录42.用抛物线插值法求解minf(x)=8x3—2x2—7x+3,已知初始单谷区间["'"]=[0,2]注=。.。。1.抛物线插值法的计算步骤：甲(t)〈甲(t),所以相对八来说t是好点，故划掉区间［t,t］,保留［t,t］为新区间,U U U乙■ JLU故置t=t，中(t)=中(t)，t=t,甲(t)=p(t)，t保持不变；2 0 2 0 0 0 1中(t)＞中(t)，所以相对，来说t是好点，故划掉区间［t,t］，保留［t,t］为新区间,0 0 1 2故置"=仍(t)*("，t°与七保持不变；程序清单抛物线插值法程序见附录5实验结果运行结果为：项目三常用无约束最优化方法(一)［实验目的］编写最速下降法、Newton法(修正Newton法)的程序。［实验准备］.掌握最速下降法的思想及迭代步骤。.掌握Newton法的思想及迭代步骤；.掌握修正Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用最速下降法求minf(X)=x2+25x;,X0=［2,2］r,£=0.01最速下降法计算步骤取初始点X(0),容许误差(精度)8>0,令k=0计算P(k)=—Vf(X(k))⑶检验||p(k)||<8?若是迭代终止,取X*=X(k),否则转4。(4)求最优步长孔:minf(X(k)+Ap(k))=f(X(k)+人/(k))(一维搜索)令X(k+1)=X(k)+人p(k),令k:=k+1,转20k最速下降法的计算框图程序清单最速下降法程序见附录6实验结果运行结果为：.用Newton法求minf(X)二60一10x-4x+X2+%2一%%初始点Xo=［0,0］『，卜0.01.Newton法的计算步骤给定初始点x(0)，及精度e>°，令k=0;Vf(X(k))<e工 —…、若 ，停止，极小点为%(◎，否则转步骤(3)；FV2f(X(k))［一1 &，(k)=一［H(X(Q)］-1Vf(X(Q)i计算 ，令 ；令x(k+1)=x(k)+s(k)k=k+1转步骤(2)程序清单Newton法程序见附录7

.97656>i]j.fli(h<r<N>"».fl0eE5PressAnykeytocentinuc.用修正Newton求minf(X)=4(x+1)2+2(x-1)2+x+x+10

1 2 1 2 ,初始点X0=[0,0]『，八0.01.修正Newton的计算步骤给定初始点x(0)，及精度'>0，令k=0;若*(X('))'&，停止，极小点为x(k)，否则转步骤(3)；「V2f(X(k))I-1 s(k)=-「H(X(k))T1Vf(X(k))计算「 」，令「 」 ；用一维搜索法求以，使得f(X(k)+a(k)S(k))=minf(X(k)+aS(k))，令O>0转步骤(2)。X(k+1)=X(k)+a(k)S(k) k转步骤(2)。， ，程序清单修正Newton程序见附录8实验结果运行结果为：实验结果运行结果为：项目四常用无约束最优化方法(二)实验目的编写共轭梯度法、变尺度法(DFP法和BFGS法)程序。实验准备1.掌握共轭方向法的思路及迭代过程；2.掌握共轭梯度法的思想及迭代步骤；3.掌握DFP法和BFGS法的思想及迭代步骤。实验内容及步骤编程解决以下问题：用共轭梯度法求得min(气2+留，取初始点X广［1，1七八0.°1.共轭梯度法算法步骤给定初始点x(0)，及精度e>0；若||巧(x(0))||V•，停止，极小值点为x(0)，否则转步骤(3)；取P(0)=—W(x(0))，且置k=0；用一维搜索法求t,使得f(x(k)+1P(k))=min/G(k)+0k)),令，k k t>0x(k+1)=x(k)+1p(k)，转步骤5；

k若|阿(x(5)1g,停止，极小值点为x(k+1)，否则转步骤⑹；若k+1=n，令x(0)=x(n)，转步骤(3),否则转步骤⑺；令P(k+1)=—W(x(k+1))+人P(k)，人= J-，置k=k+1,转步骤(4)。kk |W(x(k))||2程序清单共轭梯度法程序见附录9x±=-0.0GG002目x2=-0-003300,.pl=3_QQ00@4lj.p2=6^000088[a.h]=[-4.500O00.1.5000BB3pl=S.606604.p2=0.000000.t=0_49423Sxl=-G^GGGBBfl,x2=0-000333_pl=0_000006^)^2=-B*000005dm=[-4.500000.1.503390]=0.SS66QS^p2=-0-000005^t=0_32QQ1Q:供解为咋=-0.008000,x2=-0.000000癸■的函数值为0-000000Pi'essanykeytocentinueH.用共轭梯度法求硕5)=2""-检，自定初始点，八0.01程序清单共轭梯度法程序见附录9实验结果运行结果为：xl=B_MEIEIHHB,m2=fi.HBHBI3B,pl=-0.Bl3mEl0^p2=-0_ElBH011[a,bJ=[-4.£00000,1-500000]pl--0.660080p2~-0・00001l_,t"0.329S60量忻解为xl-0.9BS@ee,x2-0.@60009^-'W.Th0.660009Prei；t：hiiljkeyLuuuiiLlnue..用DFP法求minf(X)4(\5)+(%6),初始点X°=[8,加八0-01DFP法的具体迭代步骤如下：(1)给定初始点就顷)，迭代精度气维数n

⑵置0—k，单位矩阵I一顼七计算叫痢以计算搜索方向a进行一维搜索求：:,使必％职晔)=呼』对十/kj得迭代新点竺3气舟"T检验是否满足迭代终止条件HEQlWIWE？若满足，停止迭代，输出最优解野T)-JU：；否则进行下一步。检查迭代次数，若k=n，则置必F：”，转向步骤(2)；若k<n，则进行下一步。⑺计算：矶X矶X如U)T…一'q叫二‘产土、*•、二*甲二# ，/；，二'1]'.'.e'对蜀蜀泌｝+逐)孑WO然后，置k+1-k，转向步骤(3)。DFP法的计算框图程序清单DFP法程序见附录10实验结果运行结果为：<M1Fx2>=4<Xl-5>A2t<X2-6>"2.的极』点为:f=0弟.得机心或为自变量那伯.为，xl=5x2=6Pees注enukmutncntimie_项目五常用约束最优化方法［实验目的］编写外点罚函数法、外点罚函数法的程序。［实验准备］掌握外点罚函数法的思想及迭代步骤；掌握内点罚函数法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用外点罚函数法编程计算

minf(X)=-尤+尤,Jg(X)=ln尤>0,

.[.(X)=气+%-1=0，精度e=10-5外点法的计算步骤：给定初始点X(。)，初始罚因子人⑴，放大系数c>1；允许误差e>0,设置k=1;以x(k-i)作为搜索初始点，求解无约束规划问题minf⑴+*P(x)，令X(k为所求极小点。⑶当入P(X(k))<e，则停止计算，得到点X(k)；否则，令E)=C*)，返回⑵执行。程序清单外点罚函数法程序见附录11实验结果实验结果为运行结果为：屋[x&eAMicrosoftVisualStudio\MyProjert5\zuiyouXDebugXziyuan.exef用外点法求解：minf<X>=-Xl+X2yl<X>-lnX2>-05.t.g2<X>=X1+X2-1=0结果是：x*=7.0567e-0171f<x*>=lPressanykeytocontinue用内点罚函数法编程计算min!(x+1)3+x3i2Jx-1>0,

X>0.初始点取为X0=[3,4]「，初始障碍因子取u1=10，缩小系数取c=0,1内点罚步骤：给定初始内点工(0)eS，允许误差e>0,障碍参数丫⑴，缩小系数be(0,1)，置k=1；以x(k-1)为初始点，求解下列规划问题：minf(x)+y(k)B(x)c ，令x(k)为所求极小点s.t.xeS如果y(k)B(x(k))<e，则停止计算，得到结果x(k)，否则令y(k+1)=by(k)，置k=k+1，返回(2)。内点罚计算框图程序清单内点罚函数法程序见附录12实验结果运行结果为：附录1#include<iostream.h>#include<LIMITS.H>#defineMAX20480doublefun(doublex)(returnx*x*x-2*x+1;}doubleMax_Num(doublea,doubleb)(if(a>b)returna;elsereturnb;}doubleMin_Num(doublea,doubleb)(if(a<b)returna;elsereturnb;}voidStepAdding_Search(double&a,double&b)(doublet[MAX]={0};doubleh[MAX]={0};doublef[MAX]={0};doubleresult=0;doublep=2;t[0]=0;h[0]=1;f[0]=fun(t[0]);for(intk=0;k<MAX-1;k++){t[k+1]=t[k]+h[k];f[k+1]=fun(t[k+1]);if(f[k+1]<f[k])(h[k+1]=p*h[k];result=t[k];t[k]=t[k+1];f[k]=f[k+1];}elseif(k==0)(h[k]=-h[k];result=t[k+1];}else(a=Min_Num(result,t[k+1]);b=Max_Num(result,t[k+1]);}}}intmain()(doublea=0.0;doubleb=0.0;StepAdding_Search(a,b);cout<<"该问题的根的搜索空间是：[-<<a<<","<<b<<"]\n";return0;}附录2#include<iostream.h>#include<math.h>#defineeps0.001doublefun(doublex)(returnx*x+2*x;}doublediff_fun(doublex)(return2*x+2;}doubledichotomy(doublea,doubleb)(doublec=0.0;if((diff_fun(a)<0)&&(diff_fun(b)>0))(while(true)(c=(a+b)/2;if(diff_fun(c)<0)(a=c;if(fabs((a-b))<eps)(return(a+b)/2;}}elseif(diff_fun(c)==0)(returnc;}else(b=c;if(fabs((a-b))<eps)(return(a+b)/2;}}}}}intmain()(doublea=-3.0;doubleb=5.0;doubleresult=dichotomy(a,b);cout<<"求解的结果是:"<<result<<endl;return0;}附录3#include<iostream.h>#include<math.h>#defineeps0.01doublefun(doublex)(returnx*x*x-2*x+1;}doublediff_fun(doublex)(return3*x*x-2;}doublediff_2_fun(doublex)(return6*x;}doublenewton(doublea,doubleb)(doubleresult=(a+b)/2;doublet=0.0;while(true)(t=result-(diff_fun(result)/diff_2_fun(result));if(fabs((t-result))<eps)(returnresult;}else(result=t;}}}intmain()(doublea=0.0;doubleb=3.0;doublex=newton(a,b);cout<<"求解的结果是x="<<x<<endl;cout<<"此时f(x)="<<fun(x)<<endl;return0;}附录4#include<iostream>#include<math.h>usingnamespacestd;doublefun(doublet)(return(t*t+2*t);}voidGoldensection(double(*pfun)(doublet))(intmaxflag=1000,k=1;doublea=-3,b=5,err=0.001,t,t1,t2;do(t1=a+0.618*(b-a);t2=b-0.618*(b-a);if(t1=t2){a=t2;b=t1;}else{if(t1<t2){a=t2;}else{b=t1;}}k++;}while(fabs(a-b)>err&&k<maxflag);if(k>=maxflag)cout<<endl<<"黄金分割法迭代失败!迭代次数为k="<<k<<endl;else{t=(a+b)/2;cout<<endl<<"黄金分割法迭代成功!迭代次数为k="<<k-1<<endl;cout<<"迭代结果：近似根为root="<<t<<endl;cout<<"函数值近似为：f(root)="<<fun(t)<<endl;}}intmain(){Goldensection(fun);return0;}附录5#include<iostream>#include<math.h>usingnamespacestd;doublefun(doublex)return(8*x*x*x-2*x*x-7*x+3);}doublemax(doublea,doubleb)(if(a>b)returna;elsereturnb;}doublemin(doublea,doubleb)(if(a>b)returnb;elsereturna;}voidParainterpolation(double(*pfun)(doublex))(doublea=0,b=2,err=0.001,x=0,x0=1,f,f0;do(x=((x0*x0-b*b)*fun(a)+(b*b-a*a)*fun(x0)+(a*a-x0*x0)*fun(b))/(2*((x0-b)*fun(a)+(b-a)*fun(x0)+(a-x0)*fun(b)));f0=fun(x0);f=fun(x);if(f=f0){a=min(x,x0);b=max(x,x0);x0=(a+b)/2;}else{if((fun(x)-f0)*(x-x0)>0){b=max(x,x0);x0=min(x,x0);}else{a=min(x,x0);x0=max(x,x0);}}}while(fabs(x-x0)>err);x=(x+x0)/2;cout<<"迭代结果："<<endl;cout<<"近似根:"<<x<<endl;cout<<"函数值近似为:"<<fun(x)<<endl;}intmain(){Parainterpolation(fun);return0;}附录6#include<iostream.h>#include<math.h>doublelamda(doublex[2],doublep[2],doublea[2]){doublelam1,lam2;lam1=(pow(a[0],3)*x[0]*x[0]+pow(a[1],3)*x[1]*x[1]);lam2=-(pow(a[0]*x[0],2)+pow(a[1]*x[1],2));doubles;s=-lam2/(2*lam1);returns;}voidmain()cout<<"最速下降法求解最优解程序运行结果"<<endl;cout<<endl;doublelamd,x[3],a[6];doublep[2],g[2],e,y,m,n;inti=0;cout<<"请输入精度e"<<endl;cin>>e;cout<<"请输入初始点x[0],x[1]的值：\n"<<endl;cin>>m;cin>>n;x[0]=m;x[1]=n;cout<<"函数通式为f(x)=a[0]x1*x1+a[1]x2*x2+a[2]x1*x2+a[3]x1+a[4]x2+a[5]"<<endl;cout<<"请依次输入函数的系数：a[0]、a[1]、a[2]、a[3]、a[4]、a[5]:"<<endl;for(i=0;i<6;i++)cin>>a[i];p[0]=(2*a[0]*x[0]+a[2]*x[1]+a[3]);p[1]=(2*a[1]*x[1]+a[2]*x[0]+a[4]);g[0]=-p[0];g[1]=-p[1];i=0;cout<<endl;while(sqrt(g[0]*g[0]+g[1]*g[1])>e&&i<=200){lamd=lamda(x,g,a);x[0]=x[0]+lamd*g[0];x[1]=x[1]+lamd*g[1];p[0]=2*a[0]*x[0];p[1]=2*a[1]*x[1];g[0]=-p[0];g[1]=-p[1];i++;cout<<"'******************************************'"<<endl;cout<<"第"<<i<<"次迭代结果："<<endl;cout<<"p的模为："<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<"x的值"<<x[0]<<""<<x[1]<<endl;cout<<"'******************************************'"<<endl;cout<<endl;}y=(a[0]*x[0]*x[0]+a[1]*x[0]*x[1]+a[2]*x[0]*x[1]+a[3]*x[0]+a[4]*x[1]+a[5]);cout<<"此时满足精度要求的p的模为:"<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<endl;cout<<"满足精度的最优近似结果x[1],x[2]分别为："<<endl;cout<<"x[1]="<<x[0]<<endl;cout<<"x[2]="<<x[1]<<endl;cout<<endl;cout<<"满足进度要求所得的最优值为："<<endl;cout<<"minf(x)="<<y<<endl;}附录7#include<iostream.h>#include"matrix.h"//矩阵的有关运算intconstn1=2;doublef(double*x);//目标函数；voidg(double*x,double*y);//目标函数的梯度；voidmain(){inti;doublex0[n1],x1[n1],g0[n1],g1[n1],H[n1][n1],h[n1][n1],p[n1],f0,f1,temp,e;x0[0]=0;x0[1]=0;e=0.01;H[0][0]=2;H[0][1]=-1;H[1][0]=-1;H[1][1]=2;f0=f(x0);g(x0,g0);matrix_inv(H,h);matrix_mul(h,g0,p);//matrix_minus(x0,p,x1);for(i=0;i<n1;i++)x1[i]=x0[i]-p[i];g(x1,g1);f1=f(x1);temp=g1[0]*g1[0]+g1[1]*g1[1];while(temp>e)(for(i=0;i<n1;i++)(x0[i]=x1[i];g0[i]=g1[i];}f0=f1;matrix_mul(h,g0,p);for(i=0;i<n1;i++)x1[i]=x0[i]-p[i];g(x1,g1);f1=f(x1);temp=g1[0]*g1[0]+g1[1]*g1[1];}cout<<"当X取X("<<x1[0]<<","<<x1[1]<<")时,Min(f(X)="<<f(x1)<<endl;}doublef(double*x)(return60-10*x[0]-4*x[1]+x[0]*x[0]+x[1]*x[1]-x[0]*x[1];}voidg(double*x,double*y)(y[0]=2*x[0]-x[1]-10;y[1]=2*x[1]-x[0]-4;}附录8#include<iostream.h>#include<math.h>#include<iomanip.h>#definen2//正定二次函数的自变量个数doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值(inti,j;doublef=0;for(i=0;i<n;i++)//计算XA2部分(f+=pow(x[i],2)*f_xs[i];}for(;i<2*n;i++)//计算X部分(f+=x[i%n]*f_xs[i];}f+=f_xs[i];//计算常数项部分for(i=0;i<n;i++)//计算XiXj部分(for(j=i+1;j<n;j++)(f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];}}returnf;}voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1])(inti,j;for(i=0;i<n;i++)(Q[i][i]=2*f_xs[i];}for(i=0;i<n;i++)(Q[i][n]=f_xs[n+i];}for(i=0;i<n;i++)(for(j=i+1;j<n;j++)(Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];}}}voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n])//自变量初值，正定二次函数的各项系数，返回梯度的初值(inti,j;for(i=0;i<n;i++)(g0[i]=0;//清零for(j=0;j<n;j++)(if(i==j)g0[i]+=Q[i][j]*x[j];elseg0[i]+=Q[i][j]*x[j];}}for(i=0;i<n;i++)(g0[i]+=Q[i][n];}cout<<endl<<"梯度g0的值="<<endl;for(i=0;i<n;i++)cout<<setw(10)<<g0[i]<<endl;cout<<endl;}voidHesse_fun(doubleQ[n][n+1],doubleHesse[n][n])(inti,j;for(i=0;i<n;i++)(for(j=0;j<n;j++)(Hesse[i][j]=Q[i][j];}}cout<<"Hesse矩阵:"<<endl;for(i=0;i<n;i++)(for(j=0;j<n;j++)(cout<<setw(10)<<Hesse[i][j]<<"";}cout<<endl;}}intH(doubleg0[n],doublec)//判别准则：返回1结束，返回0继续迭代(doubles=0;for(inti=0;i<n;i++)(s+=pow(g0[i],2);}if(sqrt(s)<c)return1;elsereturn0;}voidinv(doublea[n][n],doublec[n][n])//求Hesse矩阵的逆矩阵(doublem[n];//辅助乘数doubleT;//存储换行临时变量doubletemp=0;doublet;//最大列主元inttap1=0,tap2=0;//最大列主元下标for(inti=0;i<n;i++)(for(intj=0;j<n;j++)(if(i==j)(c[i][j]=1;}else(c[i][j]=0;}}}for(ints=0;s<n;s++)(t=a[s][s];//赋初值for(intp=s;p<n;p++)//选取最大列主元(if(fabs(a[p][s])>t)(t=a[p][s];tap1=p;tap2=s;}}if(t==0){cout<<"列主元为零，失败！"<<endl;break;}if(tap1!=tap2)//换行(for(inti=0;i<n;i++)//s(T=a[tap1][i];a[tap1][i]=a[s][i];a[s][i]=T;T=c[tap1][i];//逆矩阵换行c[tap1][i]=c[s][i];c[s][i]=T;}}tap1=tap2=0;//置零for(intj=0;j<n;j++)//单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响)(a[s][j]=a[s][j]/t;c[s][j]=c[s][j]/t;}for(inti=0;i<n;i++)//消元(if(i!=s)(m[i]=a[i][s]/a[s][s];for(intj=0;j<n;j++)//注意：逆矩阵每列都消元，而原矩阵则不受影响(a[i][j]=a[i][j]-m[i]*a[s][j];c[i][j]=c[i][j]-m[i]*c[s][j];}}}}cout<<endl<<endl;cout<<"Hesse矩阵的逆矩阵为:"<<endl;for(i=0;i<n;i++)(for(intj=0;j<n;j++)(cout<<setw(10)<<c[i][j]<<"";}cout<<endl;}}voidxiang_cheng(doublea[n][n],doubleb[n],doublec[n])//n*n矩阵乘以n*1矩阵(inti,j,k;if(n==n)(cout<<"两矩阵阶数相符可以相乘！"<<endl;for(i=0;i<n;i++)(c[i]=0;}for(k=0;k<n;k++)//c矩阵的第k行(for(j=0;j<n;j++)(c[k]+=a[k][j]*b[j];}}cout<<"Hesse矩阵的逆矩阵和梯度向量相乘，两矩阵相乘结果为:"<<endl;for(i=0;i<n;i++)(cout<<setw(10)<<c[i]<<endl;}}elsecout<<"两个矩阵阶数不符，不能相乘"<<endl;}//开始计算minf(Xk+tPk)时的步长t的值，由于这是n元二次函数所以minf(t)是关于t>0的二次函数，先求二次方程a,b,c系数，用一阶导为零。voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3])//t[3]中返回的是a,b,c的系数值(inti,j;t[0]=t[1]=t[2]=0;t[0]=fun(p,f_xs)-f_xs[2*n];for(i=n;i<2*n;i++)(t[0]-=f_xs[i]*p[i%n];}for(i=0;i<n;i++)(t[1]+=2*f_xs[i]*x[i]*p[i];}for(;i<2*n;i++)(t[1]+=f_xs[i]*p[i%n];}for(i=0;i<n;i++)(for(j=i+1;j<n;j++)(t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);}}t[2]=fun(x,f_xs);cout<<endl<<"二次函数minf(Xk)的二次项系数："<<t[0]<<"—次项系数："<<t[1]<<"常数项系数："<<t[2]<<endl<<endl;}//二次函数一阶导为零计算t的值，t>0doublet_c(doublet[3])(cout<<"用一阶导求得步长为：t="<<-t[1]/t[0]/2<<endl;return-t[1]/t[0]/2;}voidmain()(doublef_xs[n+n+1+(n-1)*n/2]={4,2,9,-3,16,0};//正定二次函数的各项系数，第一个为：XS2系数，第二个为：X2A2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；//n元的系数存放类推doublex[n]={0,0};//函数自变量初值doublef0;//函数值doubleg0[n];//梯度的值doubleQ[n][n+1];//正定二次函数对应的实对称矩阵doublec=0.01;//H准则限值doubleHesse[n][n];//Hesse矩阵doubleinv_Hesse[n][n];//Hesse矩阵的逆doublet[3];//返回求minf()时t的二次函数的a,b,c的系数值doublet_bc;//步长doublep[n];//保存矩阵相乘结果_加负号后变成下降方向inti,tap=0;//迭代次数Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵ok!f0=fun(x,f_xs);//求函数初值ok!D_fun(x,Q,g0);//返回梯度的初值while(H(g0,c)==0)(Hesse_fun(Q,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变inv(Hesse,inv_Hesse);//返回Hesse矩阵的逆矩阵xiang_cheng(inv_Hesse,g0,p);for(i=0;i<n;i++)p[i]*=(-1);abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值，t_bc=t_c(t);//求一阶导来计算tfor(i=0;i<n;i++)(x[i]=x[i]+t_bc*p[i];}tap++;f0=fun(x,f_xs);D_fun(x,Q,g0);}cout<<"修正newton法"<<endl;cout<<"函数f(x1,x2)=4(X1+1)A2+2(X2-1)A2+X1+X2+10.的极小点为：’'<<"f=''<<f0<<endl;cout<<"自变量取值为："<<endl;for(i=0;i<n;i++)(cout<<"x"<<i+1<<"="<<x[i]<<endl;}cout<<"迭代次数为："<<tap<<endl;}附录9#include<stdio.h>#include<math.h>#defineN10#defineepspow(10,-6)doublef(doublex[],doublep[],doublet)(doubles;s=pow(x[0]+t*p[0],2)+4*pow(x[1]+t*p[1],2);returns;}voidsb(double*a,double*b,doublex[],doublep[])(doublet0,t1,t,h,alpha,f0,f1;intk=0;t0=2.5;/*初始值*/h=1;/*初始步长*/alpha=2;/*加步系数*/f0=f(x,p,t0);t1=t0+h;f1=f(x,p,t1);while(1)( if(f1<f0)(h=alpha*h;t=t0;t0=t1;f0=f1;k++; }else( if(k==0) (h=-h;t=t1;}else ( *a=t<t1?t:t1;*b=t>t1?t:t1;break;}}t1=t0+h;f1=f(x,p,t1);}}doublehjfg(doublex[],doublep[])(doublebeta,t1,t2,t;doublef1,f2;doublea=0,b=0;double*c,*d;c=&a,d=&b;sb(c,d,x,p);/*调用进退法搜索区间*/printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf'',x[0],x[1],p[0],p[1]);printf("\n[a,b]=[%lf,%lf]'',a,b);beta=(sqrt(5)-1.0)/2;t2=a+beta*(b-a);f2=f(x,p,t2);t1=a+b-t2;f1=f(x,p,t1);while(1)(if(fabs(t1-t2)<eps)break;else( if(f1<f2)( t=(t1+t2)/2;b=t2;t2=t1;f2=f1;t1=a+b-t2;f1=f(x,p,t1);}else( a=t1;t1=t2;f1=f2;t2=a+beta*(b-a);f2=f(x,p,t2);}}}t=(t1+t2)/2;returnt;}voidgtd()(doublex[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0;inti,k,n;printf("请输入函数的元数值n=");scanf("%d",&n);printf("\n请输入初始值:\n");for(i=0;i<n;i++)scanf("%lf'',&x[i]);f0=f(x,g,t);g[0]=2*x[0];g[1]=50*x[1];mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/if(mod1>eps)(p[0]=-g[0];p[1]=-g[1];k=0;while(1)(t=hjfg(x,p);/*调用黄金分割法求t的值*/printf("\np1=%lf,p2=%lf,t=%lf'',p[0],p[1],t);x[0]=x[0]+t*p[0];x[1]=x[1]+t*p[1];g[0]=2*x[0];g[1]=50*x[1];/*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf",x[0],x[1],g[0],g[1]);*/mod2=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/if(mod2<=eps)break;else (if(k+1==n) (g[0]=2*x[0];g[1]=50*x[1];p[0]=-g[0];p[1]=-g[1];k=0;}else(nanda=pow(mod2,2)/pow(mod1,2);printf("\nnanda=%lf,mod=%lf'',nanda,mod2);p[0]=-g[0]+nanda*p[0];p[1]=-g[1]+nanda*p[1];mod1=mod2;k++;}}printf("\n ");}}printf("\n最优解为x1=%lf,x2=%lf",x[0],x[1]);printf("\n最终的函数值为％lf",f(x,g,t));}main()(gtd();}附录10#include<iostream.h>#include<math.h>#include<iomanip.h>intconstn=2;//正定二次函数的自变量个数doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2]);//输入变量为函数自变量初值voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1]);voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n]);intH(doubleg0[n],doublec);//判别准则：返回1结束，返回0继续迭代voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3]);doublet_c(doublet[3]);//二次函数一阶导为零计算t的值，t>0voidmain()(doublef_xs[n+n+1+(n-1)*n/2]={4,1,-40,-12,136,0};doublex[n]={8,9};//函数自变量初值doublef0;//函数值doubleg0[n];//梯度的值doubleQ[n][n+1];//求梯度处值设置的中间变量，包含两部分：正定二次函数对应的实对称矩阵，还有一次项系数doublec=0.01;//H准则限值doublet[3];//返回求minf()时t的二次函数的a,b,c的系数值doublet_bc;//步长doublep[n];//保存下降方向doubleH0[n][n];//保存模拟Hesse矩阵的逆doubley[n];//y(k)=g0(k+1)-g0(k)doubles[n];//s(k)=X(k+1)-X(k)doubles_temp[n][n]={0};//计算保存矩阵doubles_temp2[n][n]={0};doubles_temp3[n][n]={0};doubles_tl[n]={0};doubletemp;//临时值inti,j,k,flag=0,tap=0;//迭代次数Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵f0=fun(x,f_xs);//求函数初值D_fun(x,Q,g0);//返回梯度的初值do{for(i=0;i<n;i++)//给H0[n][n]的处值赋单位矩阵{for(j=0;j<n;j++){if(i==j)H0[i][j]=1;elseH0[i][j]=0;}}for(i=0;i<n;i++)p[i]=(-1)*g0[i];k=0;//step2;do{abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值，t_bc=t_c(t);//求一阶导来计算tfor(i=0;i<n;i++){x[i]=x[i]+t_bc*p[i];s[i]=t_bc*p[i];//保存计算之值X(k+1)-X(k)}for(i=0;i<n;i++)y[i]=g0[i];//保存之类的f0=fun(x,f_xs);D_fun(x,Q,g0);//step3;for(i=0;i<n;i++)y[i]=g0[i]-y[i];//保存计算g0(k+1)-g0(k)if(H(g0,c)==0)//即不满足小于c(if(k!=n)(//y//shavedone!temp=0;//初值for(i=0;i<n;i++)temp+=s[i]*y[i];for(i=0;i<n;i++)(for(j=0;j<n;j++)(s_temp[i][j]=s[i]*s[j]/temp;}}//求出S(k)S(k)t/S(k)tykfor(i=0;i<n;i++)//初值s_tl[i]=0;for(i=0;i<n;i++)(for(j=0;j<n;j++)(s_tl[i]+=y[j]*H0[j][i];}}temp=0;//初值for(i=0;i<n;i++)(temp+=s_tl[i]*y[i];}//这时s_tl[n]和s_temp2[n][n]可以利用for(i=0;i<n;i++)//初值s_tl[i]=0;for(i=0;i<n;i++)(for(j=0;j<n;j++)(s_tl[i]+=H0[i][j]*y[j];}}for(i=0;i<n;i++)//初值for(j=0;j<n;j++)s_temp2[i][j]=0;for(i=0;i<n;i++)(for(j=0;j<n;j++)(s_temp2[i][j]=s_tl[i]*y[j];}}for(i=0;i<n;i++)//初值for(j=0;j<n;j++)s_temp3[i][j]=0;for(i=0;i<n;i++)(for(j=0;j<n;j++)for(intii=0;ii<n;ii++)(s_temp3[i][j]+=s_temp2[i][ii]*H0[ii][j];}}}for(i=0;i<n;i++)(for(j=0;j<n;j++)(H0[i][j]=H0[i][j]+s_temp[i][j]-s_temp3[i][j]/temp;}}for(i=0;i<n;i++)(temp=0;for(j=0;j<n;j++)(temp+=H0[i][j]*g0[j];}p[i]=(-1)*temp;}k++;}else(f0=fun(x,f_xs);break;}}else(flag=1;break;}tap++;}while(H(g0,c)==0);if(flag==1)//符合梯度准则跳出(flag=0;break;}}while(k==n);cout<<"DFP法"<<endl;cout<<"函数f(x1,x2)=4(X1-5)A2+(X2-6)A2.的极小点为："<<"f="<<f0<<endl;cout<<"取得极小点为自变量取值为："<<endl;for(i=0;i<n;i++)(cout<<"x"<<i+1<<"="<<x[i]<<endl;}}doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值(inti,j;doublef=0;for(i=0;i<n;i++)//计算XA2部分f+=pow(x[i],2)*f_xs[i];}for(;i<2*n;i++)//计算X部分(f+=x[i%n]*f_xs[i];}f+=f_xs[i];//计算常数项部分for(i=0;i<n;i++)//计算XiXj部分(for(j=i+1;j<n;j++)(f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];}}returnf;}voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1])//Q[n][n]是二次函数的正定矩阵，但Q的第n+1列存储一次项的系数(inti,j;for(i=0;i<n;i++)(Q[i][i]=2*f_xs[i];}for(i=0;i<n;i++)(Q[i][n]=f_xs[n+i];}for(i=0;i<n;i++)(for(j=i+1;j<n;j++)(Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];}}}voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n])//自变量初值，正定二次函数的各项系数，返回梯度的初值(inti,j;for(i=0;i<n;i++)(g0[i]=0;//清零for(j=0;j<n;j++)(if(i==j)g0[i]+=Q[i][j]*x[j];elseg0[i]+=Q[i][j]*x[j];}}for(i=0;i<n;i++)(g0[i]+=Q[i][n];}cout<<endl<<"梯度g0的值="<<endl;for(i=0;i<n;i++)cout<<setw(10)<<g0[i]<<endl;cout<<endl;}intH(doubleg0[n],doublec)//判别准则：返回1结束，返回0继续迭代(doubles=0;for(inti=0;i<n;i++)(s+=pow(g0[i],2);}if(sqrt(s)<c)return1;elsereturn0;}voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3])//t[3]中返回的是a,b,c的系数值(inti,j;t[0]=t[1]=t[2]=0;t[0]=fun(p,f_xs)-f_xs[2*n];for(i=n;i<2*n;i++)(t[0]-=f_xs[i]*p[i%n];}for(i=0;i<n;i++)(t[1]+=2*f_xs[i]*x[i]*p[i];}for(;i<2*n;i++)(t[1]+=f_xs[i]*p[i%n];}for(i=0;i<n;i++)(for(j=i+1;j<n;j++)(t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);}}t[2]=fun(x,f_xs);cout<<endl<<"二次函数minf(Xk)的二次项系数："<<t[0]<<"一次项系数："<<t[1]<<"常数项系数："<<t[2]<<endl<<endl;}doublet_c(doublet[3])(cout<<"用一阶导求得步长为：t="<<-t[1]/t[0]/2<<endl;return-t[1]/t[0]/2;}附录11#include<iostream.h>#include<math.h>intconstn=2;staticdoubleMk=1;doublefun_original(doublex[n]);//既是主函数，也是F(X,Mk)之一doubleg_fun(doublex[n]);//约束函数doubleh_fun(doublex[n]);;//等式函数doublea_e_f(doublex[n]);doublefun(doublex[n]);//让单纯形法求解此函数最小值voidpare(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3]);//在调用之前LGH先给x[n+1]赋值voidzx(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3],doublexn_1[n],doublexn_2[n]);doubleH(doublex[n+1][n],doubleLGH[3][n],doublec);//判别准则voiddanchun(doubleoriginal_point[n]);voidmain()(intk=1;inti;doubler=10;doublec=0.000001;doublex[n]={0,0.5};//给定的初始点在某个增广目标函数的求解域之中do{danchun(x);//ln(x2)>=0时F(X,Mk)=-X1+X2;无极小点，不用管这个增广目标函数if(Mk*a_e_f(x)<c)break;Mk=r*Mk;}while(1);cout<<"用外点法求解："<<endl;cout<<"minf(X)=-X1+X2"<<endl;cout<<" g1(X)=lnX2>=0"<<endl<<"s.t."<<endl;;cout<<" g2(X)=X1+X2-1=0"<<endl;cout<<"结果是："<<endl;cout<<"x*="<<endl;for(i=0;i<n;i++)cout<<x[i]<<"";cout<<endl;cout<<"f(x*)="<<fun_original(x)<<endl;}doublefun_original(doublex[n])//既是主函数，也是F(X,Mk)之一{doubley;y=-x[0]+x[1];returny;}doubleg_fun(doublex[n])//约束函数{doubley;y=pow(log(x[1]),2);returny;}doubleh_fun(doublex[n])//等式函数{doubley;y=x[0]+x[1]-1;returny;}doublea_e_f(doublex[n]){doubley;y=Mk*(pow(h_fun(x),2)+pow(g_fun(x),2));rsurny」doub-efun(doub-exmxAh-^^患游混稿兵因肄蜘AMW{doub-ey-ynfun—0r◎na-(X)+a—e—f(X)-rsurny」)voidpare(doub-exm+l=nLdoub-eLGHSS-doub-eLGH」aps)/7m^^N登LGH沛^兰旱旦舜园{mr-h-■doub-ewn+二for(H'o_.An+l_.++)y三nfun(x三)Wy三孔多^三三三二三三三三三二三三三三三二三二2doub-eremps」rempovyo;for(H'o_.An+l_.++)wempo^-^蜘三壁wXLiempE^-^近同xST*{一f(y=Anrempo){rempovy三」rempEn—■LGH」apovremPE7AE部克苗for(H'o_.An_.++)LGHo=vx=nu~emPE)S;for(H'o_.An+l_.++)wempo^^蜘汁壁wXHiempE^-^近同xS-FM{一f(y三vnrempo){rempovy三」rempEn—■LGH」apsnremPEWfor(H'o_.An_.++)LGHs三nx=nu~emPE)S;rempovyo;mLGH」apER。){rempovyE;for(H'o_.An+1_.++^卸洋籽汁园{一f(innLGH」apEvy一f(H'nremPE)continue-一f(y三vnrempo){rempovy三」rempEn—■LGH—rapEAempE-for(H'o_.An_.++)}voidzx(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3],doublexn_1[n],doublexn_2[n])(inti,j;for(i=0;i<n;i++)xn_1[i]=xn_2[i]=0;for(i=0;i<n+1;i++)//返回重心(if(i==LGH_tap[2])//轮流加n次(去掉一次)continue;for(j=0;j<n;j++)(xn_1[j]+=x[i][j];}}for(i=0;i<n;i++)xn_1[i]/=n;for(i=0;i<n;i++)//返回反射点xn_2[i]=2*xn_1[i]-LGH[2][i];}doubleH(doublex[n+1][n],doubleLGH[3][n],doublec)//判别准则(inti;doubles=0;for(i=0;i<n+1;i++)//小心i从0到n即(i=0;i<n+1;i++)(s+=pow(fun(x[i])-fun(LGH[0]),2);//s+=pow(fun(x[i][n])-fun(LGH[0][n]),2);//////////////}if(s<=c)return1;elsereturn0;}voiddanchun(doubleoriginal_point[n])(doublex[n+1][n]={0};doublec=0.000001;doublea=2;doubleb=0.5;doubleLGH[3][n]={0};doubleLGH_tap[3]={0};doublexn_1[n]={0};doublexn_2[n]={0};doublexn_3[n]={0};doublexn_4[n]={0};inti,j,tap=0;for(i=0;i<n;i++)x[0][i]=original_point[i];for(i=1;i<n+1;i++)//假设只给出一个初始点X0(1,1),求另外两个初始点{for(j=0;j<n;j++){if(j+1==i)x[i][j]=x[0][j]+1;else}pare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它们在x[n+1][n]中的下标zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射点X(n+2)while(H(x,LGH,c)==0)(if(fun(xn_2)<fun(LGH[0]))//if(fun(xn_2)<fun(LGH[0][n]))//////////////////(for(i=0;i<n;i++)(xn_3[i]=xn_1[i]+a*(xn_2[i]-xn_1[i]);//扩X}if(fun(xn_3)<fun(LGH[0]))(for(i=0;i<n;i++)(x[int(LGH_tap[2])][i]=LGH[2][i]=xn_3[i];}}else(for(i=0;i<n;i++)(x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}}}else(if(fun(xn_2)<fun(LGH[1]))(for(i=0;i<n;i++)(x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}//gotostep2;}else(if(fun(xn_2)<fun(LGH[2]))//elseif(fun(xn_2)<fun(LGH[2][n]))(for(i=0;i<n;i++)(x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}}else//means(fun(xn_2)>=fun(LGH[2][n]))(for(i=0;i<n;i++)(xn_4[i]=xn_1[i]+b*(xn_2[i]-xn_1[i]);//压缩}if(fun(xn_4)<fun(LGH[2]))(for(i=0;i<n;i++)x[int(LGH_tap[2])][i]=LGH[2][i]=xn_4[i];}//gotostep2;}else(for(i=0;i<n+1;i++)(if(i==LGH_tap[0])continue;for(j=0;j<n;j++)(x[i][j]=(LGH[0][j]+x[i][j])/2;}}}//gotostep2}}}//That'sit!!!pare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它们在x[n+1][n]中的下标zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射点X(n+2)tap++;}for(i=0;i<n;i++)original_point[i]=LGH[0][i];//把在某个Mk值下增广目标函数的最优值返回}附录12#include<iostream.h>#include<math.h>#include<iomanip.h>constintn=2;staticdoubleRk=10;//Rk是惩罚因子staticintp=1;doublefun_original(doublex[n]);//既是主函数，也是F(X,Rk)之一doublea_e_f(doublex[n]);doublefun(doublex[n]);//让单纯形法求解此函数最小值voidg_fun(doublex[n],doubleg0[n]);//牛顿法子函数开始,求出与Rk对应的梯度向量.voidHesse_fun(doublex[n],doubleHesse[n][n]);//求出与Rk对应的Hesse矩阵intH(doubleg0[n],doublec);voidinv(doublea[n][n]