最优化方法课程实验报告

项目维搜索算法(一)［实验目的］编写加步探索法、对分法、Newton法的程序。［实验准备］1.掌握一维收搜索中搜索区间的加步探索法的思想及迭代步骤；2.掌握对分法的思想及迭代步骤；3.掌握Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用加步探索法确定一维最优化问题min9(t)=13—2t+1t>0的搜索区间，要求选取10=°，h0=1,a=2.加步探索法算法的计算步骤:⑴选取初始点t£［0，+8)(或，£U［°，t］)计算甲=9(t)给出初始步长h〉0I/J>1-1八、、00max,I牙00♦l口a*h0A0,加步系数以>1,令k=0。比较目标函数值.令t=t+h,计算9=9(t)若9<9，转(3)，否k+1kkk+1k+1k+1k贝膊专(4)。(o\加十捽步中甘Ah=ah同时At=t,t=t,k=k+1姓加大探索步长♦令k+1k，同时，令k'kk+1,,转(2)。反向探索.若k=0,转换探索方向,令气=-气,t=\，转(2)。否则，停止迭代，令a=代，令a=min{t，t}，b=max{t，t}加步探索法算法的计算框图程序清单加步探索法算法程序见附录1实验结果运行结果为：I'D:\Pnogr3mFile&(croso-HVisualStudio\k/lyPr该问题的根的搜索空间是：［^31PressanyIce-tocontinue.2.用对分法求解min中(t)=t(t+2),已知初始单谷区间［a，b］=［—3,5］,要求按精度£=。.3,£=0.001分别计算.对分法迭代的计算步骤：确定初始搜索区间［a,b］,要求中'(a)<0,中啊>0。计算［a,b］的中点c=1(a+b)2'若甲'(c)V0，则a=c,转(4);若甲'(c)=0，则t*=c,转(5);若甲'(c)>0，则b=c，转(4).

1/八(4)右|a-bl<e，则t*=-(a+b),转(5);否则转(2).^2⑸打印t*，结束程序清单对分法程序见附录2实验结果运行结果为：I'D:\PnogramFiles(x86AMicrc£oftVisualStudic\MyProject世解的结果是：TPressanykmytocontinue3.用Newton法求解min中(t)=13—2t+1/已知初始单谷区间[a，b]=[0，1]，要求精度8=0.01.Newton法的计算步骤

⑴确定初始搜索区间［。,仞，要求甲’(a)<0，平’(b)>0(2)选定t0⑶计算t='。—甲甲"(t0)若11-10l^e，则t0=t，转(3);否则转(5).打印t，中(t)，结束.Newton法的计算框图程序清单Newton法程序见附录3实验结果运行结果为:项目二一维搜索算法(二)［实验目的］编写黄金分割法、抛物线插值法的程序。［实验准备］掌握黄金分割法的思想及迭代步骤；掌握抛物线插值法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用黄金分割法求解min中(t)=t(t+2),已知初始单谷区间［a，b］=［—3,5］,要求精度8=0.001.黄金分割法迭代步骤:确定中(t)的初始搜索区间［a，b］.计算12=a+0.382(b-a)计算ti=a+0.618(b-a)若1ti-t2l<e,则打印t*=1^2，结束；否则转(5).判别是否满足甲V甲：若满足，则置a=t，t=t，甲=平1222121然后转(3);否则，置b=t，t=t，甲=甲，t=a+。(b—a)，甲=甲(t)11212222然后转(4).黄金分割法的计算框图：程序清单黄金分割法程序见附录4实验结果运行结果为：2.用抛物线插值法求解minf(x)=8x3—2x2—7x+3已知初始单谷区间［",』］-心二顷.抛物线插值法的计算步骤：甲(t)〈甲(t),所以相对t来说t是好点，故划掉区间［t,t］,保留［t,t］为新区TOC\o"1-5"\h\z000210间,故置t=t，甲(t)=甲(t),t=t，甲(t)=甲(t),t保持不变;2020001甲(t)＞甲(t),所以相对t来说t是好点，故划掉区间［t,t］,保留［t,t］为新区间，0012故置t1=t，^()=中(t1),t0与12保持不变；程序清单抛物线插值法程序见附录5实验结果运行结果为:项目三常用无约束最优化方法(一)［实验目的］编写最速下降法、Newton法(修正Newton法)的程序。［实验准备］.掌握最速下降法的思想及迭代步骤。.掌握Newton法的思想及迭代步骤;.掌握修正Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：用最速下降法求minf(X)=x2+25x;，X0=［2,2］r，£=0.01.最速下降法计算步骤(1)取初始点X(0)，容许误差(精度)£>0，令k=0(2)计算p(k)=-Vf(X(k))(3)检验||p(k<e?若是迭代终止,取X*=X(k),否则转4。(4)求最优步长孔：minf(X(k)+Xp(k))=f(X(k)+Rp(k))(一维搜索)令X(k+i)=X(k)+Xp(k),令k:=k+1,转20k最速下降法的计算框图最速下降法程序见附录6实验结果运行结果为:2.用Newton法求minf(X)二60一10x-4x+X2+%2一%%初始点X0=[0，0]T，8=0.01Newton法的计算步骤(1)给定初始点x(0),及精度8>0，令k=0；(2)若HVf(X(k))1~8,停止，极小点为％(k)，否则转步骤(3);(3)计算N2f(X(k))]T,令s(k)=-[H(X(k))]TVf(X(k))；令%(k+i)=%(k)+s(k),k=k+1,转步骤(2)。程序清单Newton法程序见附录7实验结果运行结果为:Frcssanykeytocontinuc3.用修正Newton求minf(X)=4(%+1)2+2(%一1)2+%+%+10,

1212初始点X0=[0，0]T，8=0.01.修正Newton的计算步骤(1)给定初始点%(0),及精度8>0,令k=0;(2)若HVf(X(k叫~8,停止，极小点为％(k)，否则转步骤(3);

(3)计算N^f(X(人))J',令S(人)=-[H(X(k))J1巧(X(k))；(4)用一维搜索法求a,使得f(X(k)+a(Sk沪minfX^a)Sk,)令a>0X(k+1)=X(k)+a(k)S(k),k=k+1,转步骤(2)。程序清单修正Newton程序见附录8实验结果运行结果为：T'DeprogramFiles；x8c)\Micro=oftVis-uaJStudicVMyPrcyects^zuiyciijXC^bLg^zryuan.exe'向量相乘，两矩阵相乘结果为、二状函薮FM湖农顼系数'S-iS^—农顼系数'-12.375常数项系数・16椽盅newt<m法函数ESL或〉T就』*1十2<X2-1)A2#S1■*«2*10一的极小点为i自堡童取值为:向量相乘，两矩阵相乘结果为、>rl=-lA25直代次数为：1Pi'bisean事keytocointinLLB项目四常用无约束最优化方法（二）实验目的编写共轭梯度法、变尺度法（DFP法和BFGS法）程序。实验准备1.掌握共轭方向法的思路及迭代过程；掌握共轭梯度法的思想及迭代步骤；3.掌握DFP法和BFGS法的思想及迭代步骤。实验内容及步骤编程解决以下问题：1.用共轭梯度法求得min(时+4弋),取初始点％=[1，1]r,e=0.01.共辆梯度法算法步骤(1)给定初始点x(0),及精度e>0;(2)若|阿(x(0))1<£,停止,极小值点为X(0),否则转步骤(3);(3)取p(0)=—Nf(x(0)),且置k=0；(4)用一维搜索法求t，使得f(xk)+tpk)=minfx+tp)令，X(k+1)=x(k)+tp(k),转步骤5；k(5)若||Vf(x(k+1))|<e,停止,极小值点为x(k+1),否则转步骤(6);(6)若k+1=n，令x(0)=x(n),转步骤(3),否则转步骤(7);Vf(x(k+1))12(7)令p(k+1)=—Vf(x(k+1))+人p(k),人=,置k=k+1,转步骤kk||Vf(x(k))||2。程序清单共轭梯度法程序见附录9实验结果运行结果为：

D:\ProgramFiles(x86)\IMicro5oftVisualStudio\MyPrqject5\zuiyou\Debug\ziyuan.exe"阿输入函数的元数值忡请输入初始值：11=-B.000000,pl=fl.013131304^p2=B_ra.b：l=I-4.500BB0.1-500000]=B.4?4238xl=^0.060000,x2-000000,pi=0・000090^p2=-0.000005[a^b]-[-4.566660,1.500000]pl-0-666606,p2—6.600005020913最代解为xi-g.GeGeGo^xE-e.eoeoea最教色匡苏值为白.GGGOOOrressanjpkeytocontimiE：—2.用共轭梯度法求min/（X）=2^+X2-x1七，自定初始点，e=0.01.程序清单共轭梯度法程序见附录9实验结果运行结果为:rr"D:\ProgramiFiles(K86)\MicrD5aftVisualStudioVMyPT^jertsXzui^ouXDebugXziyuan.exe8清输A函分的打数道n=W_.j—一—清输A函分的打数道n=W睛输入初始值:Y1=n一MMWMMMrY2=H.1买=-M.MWMMMMP=-rt.MMMM11Pt=rt.僵技由函数0-B300G0Fressany虹母tncontinue.3.用DFP法求minf(X)=4(%—5)2+(气-6)2,初始点X0=[8,9]…0.01.DFP法的具体迭代步骤如下:(1)给定初始点对),迭代精度，，维数n(2)置0—k,单位矩阵I—,计算矽•。(3)计算搜索方向-•(4)进行一维搜索求危，使得迭代新点济3气舟"7(5)检验是否满足迭代终止条件H'l\lII?若满足,停止迭代,输出最优解野T),_厂，；否则进行下一步。(6)检查迭代次数，若k=n，则置x^-xin),转向步骤(2);若k<n，则进行下一步。(7)计算:-必'TA'然后，置k+1-k,转向步骤（3）。DFP法的计算框图程序清单DFP法程序见附录10实验结果运行结果为:DFF法函数”心“2〉=4牌1-5〉七+牌2-6〉七.的极小点为:f=0刑得板小点为日变量卅伯.为，xl=5^2=6点目nukeutncritinu皓一项目五常用约束最优化方法［实验目的］编写外点罚函数法、外点罚函数法的程序。［实验准备］掌握外点罚函数法的思想及迭代步骤；掌握内点罚函数法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：1.用外点罚函数法编程计算minf(X)=-尤+尤，Jg(X)=ln尤>0,.［.(X)=气+%-1=0，精度e=10-5.外点法的计算步骤：(1)给定初始点x(0),初始罚因子从1),放大系数c>1;允许误差e>0,设置k=1;(2)以x(k-1)作为搜索初始点，求解无约束规划问题minf(x)+人P(x)，令X(k)为所求极小点。(3)当人P(x(们)<°,则停止计算,得到点x(k)；否则，令Xk+1)=ck(k),返回(2)执行。程序清单外点罚函数法程序见附录11实验结果实验结果为运行结果为：：s(x86hMicroso-HVisualStudio\MyProjects^zuiyou\Debug\ziyLjan.exe'担外点法求解：tiinf<X>=-X1+X2gl<X>=lnX2>=0s.t.g2<X>=Xl+X2-l=0结果是：x*=7.0567e-0171f<x*>=1Pressanykeytocontinue

2.用内点罚函数法编程计算min[3(七+1)3+%」，0.1JX-1>0,

X>0.0.1初始点取为X0=[3,4H，初始障碍因子取"1=10，缩小系数取C内点罚步骤：给定初始内点X(0)GS，允许误差e>0,障碍参数Y⑴，缩小系数bG(0,1)，置k=1；以X(k-1)为初始点，求解下列规划问题：minf(x)+y(k)B(x)人，令X(k)为所求极小点st.XGS(3)如果y(k)B(X(k))<e，则停止计算，得到结果X(k),(4)否则令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[i]=-p[i];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]分别为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++)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++)//计算X人2部分{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};//正定二次函数的各项系数,第一个为：X1人2系数,第二个为：X2人2系数,第三个为：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)人2+2(X2-1)人2+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,T2,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)人2+(X2-6)人2.的极小点为："<<"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++)//计算X人2部分{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]);//让单纯形法求解此函数最小值voidcompare(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));returny;}doublefun(doublex[n])//让单纯形法求解此函数最小值{doubley;y=fun_original(x)+a_e_f(x);returny;}voidcompare(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3])//在调用之前LGH先给x[n+1]赋值{inti;doubley[n+1];for(i=0;i<n+1;i++)y[i]=fun(x[i]);//y[i]=fun(x[i][n]);//////////////////////////////////////////////////?doubletemp[2];temp[0]=y[0];for(i=0;i<n+1;i++)//temp[0]存储最小的值XL,temp[1]存储对应x的下标{if(y[i]<=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[0]=temp[1];//记录标记for(i=0;i<n;i++)LGH[0][i]=x[int(temp[1])][i];for(i=0;i<n+1;i++)//temp[0]存储最大的值XH,temp[1]存储对应x的下标{if(y[i]>=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[2]=temp[1];//for(i=0;i<n;i++)LGH[2][i]=x[int(temp[1])][i];temp[0]=y[0];if(LGH_tap[2]==0){temp[0]=y[1];}for(i=0;i<n+1;i++)//寻找次大值{if(i==LGH_tap[2])//if(i==temp[1])continue;if(y[i]>=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[1]=temp[1];for(i=0;i<n;i++)LGH[1][i]=x[int(temp[1])][i];}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;elsex[i][j]=x[0][j];}}compare(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]);//扩张}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!!!compare(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对应的Hes