最优化方法课程实验报告.doc

项目一 一维搜索算法（一）实验目的编写加步探索法、对分法、Newton法的程序。实验准备1掌握一维收搜索中搜索区间的加步探索法的思想及迭代步骤；2掌握对分法的思想及迭代步骤；3掌握Newton法的思想及迭代步骤。实验内容及步骤 编程解决以下问题：1用加步探索法确定一维最优化问题的搜索区间，要求选取加步探索法算法的计算步骤： (1)选取初始点，计算给出初始步长，加步系数，令。 (2) 比较目标函数值令，计算 ，若，转(3)，否则转(4)。(3) 加大探索步长令，同时，令，转(2)。 (4) 反向探索若，转换探索方向，令，转（2）。否则，停止迭代，令。加步探索法算法的计算框图程序清单加步探索法算法程序见附录1实验结果运行结果为：2用对分法求解，已知初始单谷区间，要求按精度，分别计算对分法迭代的计算步骤：(1)确定初始搜索区间，要求。 (2) 计算的中点 (3) 若，则 ，转(4)；若，则，转(5)；若，则 ，转(4) (4) 若，则，转(5)；否则转(2)(5) 打印，结束对分法的计算框图程序清单对分法程序见附录2实验结果运行结果为：3用Newton法求解，已知初始单谷区间，要求精度Newton法的计算步骤(1) 确定初始搜索区间，要求 (2) 选定(3) 计算 (4) 若 ，则，转（3）；否则转（5） (5) 打印 ，结束 Newton法的计算框图 程序清单Newton法程序见附录3实验结果运行结果为：项目二 一维搜索算法（二）实验目的编写黄金分割法、抛物线插值法的程序。实验准备1掌握黄金分割法的思想及迭代步骤；2掌握抛物线插值法的思想及迭代步骤。实验内容及步骤编程解决以下问题：1用黄金分割法求解，已知初始单谷区间，要求精度黄金分割法迭代步骤： (1) 确定的初始搜索区间 (2) 计算(3) 计算(4) 若，则打印，结束；否则转(5) (5) 判别是否满足：若满足，则置 然后转(3)；否则，置然后转(4)黄金分割法的计算框图：程序清单黄金分割法程序见附录4实验结果运行结果为： 2用抛物线插值法求解，已知初始单谷区间抛物线插值法的计算步骤：（1） ，所以相对来说是好点，故划掉区间，保留为新区间，故置，保持不变；（2） ，所以相对来说是好点，故划掉区间，保留为新区间，故置，与保持不变； 程序清单 抛物线插值法程序见附录5实验结果运行结果为：项目三 常用无约束最优化方法（一）实验目的编写最速下降法、Newton法（修正Newton法）的程序。实验准备1掌握最速下降法的思想及迭代步骤。2掌握Newton法的思想及迭代步骤；3掌握修正Newton法的思想及迭代步骤。实验内容及步骤编程解决以下问题：1用最速下降法求最速下降法计算步骤（1）（2）（3）（4）最速下降法的计算框图程序清单最速下降法程序见附录6实验结果运行结果为：2用Newton法求，初始点Newton法的计算步骤（1）给定初始点，及精度，令；（2）若，停止，极小点为，否则转步骤（3）；（3）计算，令；令，转步骤（2）。程序清单Newton法程序见附录7实验结果 运行结果为：3用修正Newton求，初始点修正Newton的计算步骤（1）给定初始点，及精度，令；（2）若，停止，极小点为，否则转步骤（3）；（3）计算，令；（4）用一维搜索法求，使得，令，转步骤（2）。程序清单修正Newton程序见附录8实验结果运行结果为：项目四 常用无约束最优化方法（二）实验目的编写共轭梯度法、变尺度法（DFP法和BFGS法）程序。实验准备1掌握共轭方向法的思路及迭代过程；2掌握共轭梯度法的思想及迭代步骤；3掌握DFP法和BFGS法的思想及迭代步骤。实验内容及步骤编程解决以下问题：1 用共轭梯度法求得，取初始点，共轭梯度法算法步骤（1） 给定初始点，及精度；（2） 若，停止，极小值点为，否则转步骤（3）；（3） 取，且置；（4） 用一维搜索法求，使得，令，转步骤5；（5） 若，停止，极小值点为，否则转步骤（6）；（6） 若，令，转步骤（3），否则转步骤（7）；（7） 令，置，转步骤（4）。程序清单共轭梯度法程序见附录9实验结果运行结果为：2 用共轭梯度法求，自定初始点，程序清单共轭梯度法程序见附录9实验结果运行结果为：3用DFP法求，初始点DFP法的具体迭代步骤如下：（1）给定初始点，迭代精度，维数n（2）置0k，单位矩阵I，计算。（3）计算搜索方向（4）进行一维搜索求，使得迭代新点（5）检验是否满足迭代终止条件？若满足，停止迭代，输出最优解,；否则进行下一步。（6）检查迭代次数，若k=n，则置，转向步骤(2)；若k1；允许误差e0,设置k=1； （2）以x(k-1)作为搜索初始点，求解无约束规划问题，令x(k)为所求极小点。 （3）当，则停止计算，得到点x(k)；否则，令，返回（2）执行。程序清单外点罚函数法程序见附录11实验结果实验结果为运行结果为：2用内点罚函数法编程计算初始点取为，初始障碍因子取，缩小系数取内点罚步骤：（1） 给定初始内点，允许误差e0，障碍参数，缩小系数，置k=1；（2） 以为初始点，求解下列规划问题：，令为所求极小点（3） 如果，则停止计算，得到结果，（4） 否则令，置k=k+1，返回（2）。内点罚计算框图程序清单内点罚函数法程序见附录12实验结果运行结果为：附录1#include #include #define MAX 20480 double fun(double x) return x*x*x-2*x+1; double Max_Num(double a,double b) if(ab) return a; else return b; double Min_Num(double a,double b) if(ab) return a; else return b; void StepAdding_Search(double &a,double &b) double tMAX=0; double hMAX=0; double fMAX=0; double result=0; double p=2; t0=0; h0=1; f0=fun(t0); for(int k=0;kMAX-1;k+) tk+1=tk+hk; fk+1=fun(tk+1); if(fk+1fk) hk+1=p*hk; result=tk; tk=tk+1; fk=fk+1; else if(k = 0) hk=-hk; result=tk+1; else a=Min_Num(result,tk+1); b=Max_Num(result,tk+1); int main() double a = 0.0; double b= 0.0; StepAdding_Search(a,b); cout该问题的根的搜索空间是: a,bn; return 0; 附录2#include #include #define eps 0.001 double fun(double x) return x*x+2*x; double diff_fun(double x) return 2*x+2; double dichotomy(double a,double b) double c = 0.0; if(diff_fun(a) 0) while(true) c = (a + b)/2; if(diff_fun(c) 0) a=c; if(fabs(a-b) eps) return (a+b)/2; else if(diff_fun(c) = 0) return c; else b = c; if(fabs(a-b) eps) return (a+b)/2; int main() double a = -3.0; double b = 5.0; double result = dichotomy(a,b); cout求解的结果是:resultendl; return 0; 附录3#include #include #define eps 0.01 double fun(double x) return x*x*x-2*x+1; double diff_fun(double x) return 3*x*x-2; double diff_2_fun(double x) return 6*x; double newton(double a,double b) double result = (a+b)/2; double t = 0.0; while(true) t = result - (diff_fun(result)/diff_2_fun(result); if(fabs(t-result) eps) return result; else result =t; int main() double a = 0.0; double b = 3.0; double x= newton(a,b); cout求解的结果是 x = xendl; cout此时f(x) = fun(x)endl; return 0; 附录4#include#includeusing namespace std;double fun(double t)return (t*t+2*t);void Goldensection(double(*pfun)(double t)int maxflag=1000,k=1;double a=-3,b=5,err=0.001,t,t1,t2;dot1=a+0.618*(b-a);t2=b-0.618*(b-a);if(t1=t2)a=t2;b=t1;else if(t1err&k=maxflag)coutendl黄金分割法迭代失败!迭代次数为k=kendl; elset=(a+b)/2;coutendl 黄金分割法迭代成功!迭代次数为k=k-1endl;cout 迭代结果：近似根为root=tendl; cout 函数值近似为：f(root)=fun(t)endl; int main()Goldensection(fun);return 0;附录5#include#includeusing namespace std;double fun(double x)return (8*x*x*x-2*x*x-7*x+3);double max(double a,double b)if(ab)return a;else return b;double min(double a,double b)if(ab)return b;else return a;void Parainterpolation(double(*pfun)(double x) double a=0,b=2,err=0.001,x=0,x0=1,f,f0;dox=(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);elsea=min(x,x0);x0=max(x,x0);while(fabs(x-x0)err);x=(x+x0)/2;cout 迭代结果：endl;cout 近似根:xendl;cout 函数值近似为:fun(x)endl;int main()Parainterpolation(fun);return 0;附录6#include#includedouble lamda(double x2,double p2,double a2) double lam1,lam2; lam1=(pow(a0,3)*x0*x0+pow(a1,3)*x1*x1); lam2=-(pow(a0*x0,2)+pow(a1*x1,2); double s;s=-lam2/(2*lam1); return s; void main() cout最速下降法求解最优解程序运行结果endl; coutendl; double lamd,x3,a6; double p2,g2,e,y,m,n; int i=0; cout请输入精度ee; cout请输入初始点x0,x1的值：nm; cinn; x0=m; x1=n; cout函数通式为f(x)=a0x1*x1+a1x2*x2+a2x1*x2+a3x1+a4x2+a5endl; cout请依次输入函数的系数：a0、a1、a2、a3、a4、a5:endl; for(i=0;iai;p0=(2*a0*x0+a2*x1+a3); p1=(2*a1*x1+a2*x0+a4); g0=-p0; g1=-p1; i=0; coute&i=200) lamd=lamda(x,g,a); x0=x0+lamd*g0; x1=x1+lamd*g1;p0=2*a0*x0; p1=2*a1*x1; g0=-p0; g1=-p1; i+;cout*endl; cout第i次迭代结果：endl; coutp的模为：sqrt(g0*g0+g1*g1)endl; coutx的值x0 x1endl; cout*endl; coutendl; y=(a0*x0*x0+a1*x0*x1+a2*x0*x1+a3*x0+a4*x1+a5);cout此时满足精度要求的p的模为:sqrt(g0*g0+g1*g1)endl; coutendl;cout满足精度的最优近似结果x1,x2分别为：endl; coutx1=x0endl; coutx2=x1endl; coutendl; cout满足进度要求所得的最优值为：endl;coutminf(x)=yendl; 附录7#include#includematrix.h/矩阵的有关运算int const n1=2;double f(double *x);/目标函数;void g(double *x,double *y);/目标函数的梯度;void main()int i;double x0n1,x1n1,g0n1,g1n1,Hn1n1,hn1n1,pn1,f0,f1,temp,e;x00=0;x01=0;e=0.01;H00=2;H01=-1;H10=-1;H11=2;f0=f(x0);g(x0,g0);matrix_inv(H,h);matrix_mul(h,g0,p);/matrix_minus(x0,p,x1);for(i=0;ie)for(i=0;in1;i+)x0i=x1i;g0i=g1i;f0=f1;matrix_mul(h,g0,p);for(i=0;in1;i+)x1i=x0i-pi;g(x1,g1);f1=f(x1);temp=g10*g10+g11*g11;cout当X取X(x10,x11)时,Min(f(X)=f(x1)endl;double f(double *x)return 60-10*x0-4*x1+x0*x0+x1*x1-x0*x1;void g(double *x,double *y)y0=2*x0-x1-10;y1=2*x1-x0-4;附录8#include#include#include#define n 2/正定二次函数的自变量个数double fun(double xn,double f_xsn+n+1+(n-1)*n/2)/输入变量为函数自变量初值int i,j;double f=0;for(i=0;in;i+)/计算X2部分f+=pow(xi,2)*f_xsi;for(;i2*n;i+)/计算X部分f+=xi%n*f_xsi;f+=f_xsi;/计算常数项部分for(i=0;in;i+)/计算XiXj部分for(j=i+1;jn;j+)f+=f_xs(n+1)+n*(i+1)-i*(i+1)/2+j-i-1*xi*xj;return f;void Q_fun(double f_xsn+n+1+(n-1)*n/2,double Qnn+1)int i,j;for(i=0;in;i+)Qii=2*f_xsi;for(i=0;in;i+)Qin=f_xsn+i;for(i=0;in;i+)for(j=i+1;jn;j+)Qji=Qij=f_xs(n+1)+n*(i+1)-i*(i+1)/2+j-i-1;void D_fun(double xn,double Qnn+1,double g0n)/自变量初值，正定二次函数的各项系数，返回梯度的初值int i,j;for(i=0;in;i+)g0i=0;/清零for(j=0;jn;j+)if(i=j)g0i+=Qij*xj;else g0i+=Qij*xj;for(i=0;in;i+)g0i+=Qin;coutendl梯度g0的值=endl;for(i=0;in;i+)coutsetw(10)g0iendl;coutendl;void Hesse_fun(double Qnn+1,double Hessenn)int i,j;for(i=0;in;i+)for(j=0;jn;j+)Hesseij=Qij;coutHesse 矩阵:endl;for(i=0;in;i+)for(j=0;jn;j+)coutsetw(10)Hesseij ;coutendl;int H(double g0n,double c)/判别准则：返回1结束，返回0继续迭代double s=0;for(int i=0;in;i+)s+=pow(g0i,2);if(sqrt(s)c)return 1;else return 0;void inv(double ann,double cnn)/求Hesse矩阵的逆矩阵double mn;/辅助乘数double T;/存储换行临时变量double temp=0;double t;/最大列主元int tap1=0,tap2=0;/最大列主元下标for(int i=0;in;i+)for(int j=0;jn;j+)if(i=j)cij=1;elsecij=0;for(int s=0;sn;s+)t=ass;/赋初值for(int p=s;pt)t=aps;tap1=p;tap2=s;if(t=0)cout列主元为零，失败！endl;break;if(tap1!=tap2)/换行for(int i=0;in;i+)/sT=atap1i;atap1i=asi;asi=T;T=ctap1i;/逆矩阵换行ctap1i=csi;csi=T;tap1=tap2=0;/置零for(int j=0;jn;j+)/单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响)asj=asj/t;csj=csj/t;for(int i=0;in;i+)/消元if(i!=s)mi=ais/ass;for(int j=0;jn;j+)/注意：逆矩阵每列都消元，而原矩阵则不受影响aij=aij-mi*asj;cij=cij-mi*csj;coutendlendl;coutHesse矩阵的逆矩阵为:endl;for(i=0;in;i+)for(int j=0;jn;j+) coutsetw(10)cij ;coutendl;void xiang_cheng(double ann,double bn,double cn)/n*n矩阵乘以n*1矩阵int i,j,k;if(n=n)cout两矩阵阶数相符可以相乘！endl;for(i=0;in;i+)ci=0;for(k=0;kn;k+)/c矩阵的第k行for(j=0;jn;j+)ck+=akj*bj;coutHesse矩阵的逆矩阵和梯度向量相乘，两矩阵相乘结果为:endl;for(i=0;in;i+)coutsetw(10)ciendl;elsecout两个矩阵阶数不符，不能相乘0的二次函数，先求二次方程a,b,c系数,用一阶导为零。void abc(double xn,double pn,double f_xsn+n+1+(n-1)*n/2,double t3)/t3中返回的是a,b,c的系数值int i,j;t0=t1=t2=0;t0=fun(p,f_xs)-f_xs2*n;for(i=n;i2*n;i+)t0-=f_xsi*pi%n;for(i=0;in;i+)t1+=2*f_xsi*xi*pi;for(;i2*n;i+)t1+=f_xsi*pi%n;for(i=0;in;i+)for(j=i+1;jn;j+)t1+=f_xs(n+1)+n*(i+1)-i*(i+1)/2+j-i-1*(xi*pj+xj*pi);t2=fun(x,f_xs);coutendl二次函数minf(Xk)的二次项系数：t0 一次项系数：t1 常数项系数：t2endl0double t_c(double t3)cout用一阶导求得步长为：t=-t1/t0/2endl;return -t1/t0/2;void main()double f_xsn+n+1+(n-1)*n/2=4,2,9,-3,16,0;/正定二次函数的各项系数，第一个为：X12系数，第二个为：X22系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；/n元的系数存放类推double xn=0,0;/函数自变量初值double f0;/函数值double g0n;/梯度的值double Qnn+1;/正定二次函数对应的实对称矩阵double c=0.01;/H准则限值double Hessenn;/Hesse矩阵double inv_Hessenn;/Hesse矩阵的逆double t3;/返回求minf()时t的二次函数的a,b,c的系数值double t_bc;/步长double pn;/保存矩阵相乘结果_加负号后变成下降方向int i,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;in;i+)pi*=(-1); abc(x,p,f_xs,t);/开始计算minf(Xk+tPk)时的步长t的值，t_bc=t_c(t);/求一阶导来计算tfor(i=0;in;i+)xi=xi+t_bc*pi;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=f0endl;cout自变量取值为：endl;for(i=0;in;i+)coutxi+1=xiendl;cout迭代次数为：tapendl;附录9#include #include #define N 10 #define eps pow(10,-6)double f(double x,double p,double t) double s; s=pow(x0+t*p0,2)+4*pow(x1+t*p1,2); return s; void sb(double *a,double *b,double x,double p)double t0,t1,t,h,alpha,f0,f1; int k=0; t0=2.5; /*初始值*/ h=1; /*初始步长*/ alpha=2; /*加步系数*/ f0=f(x,p,t0); t1=t0+h; f1=f(x,p,t1);while(1) if(f1f0) h=alpha*h; t=t0; t0=t1; f0=f1; k+; else if(k=0) h=-h;t=t1; else *a=tt1?t:t1; break; t1=t0+h; f1=f(x,p,t1); double hjfg(double x,double p)double beta,t1,t2,t; double f1,f2; double a=0,b=0; double *c,*d; c=&a,d=&b; sb(c,d,x,p);/*调用进退法搜索区间*/ printf(nx1=%lf,x2=%lf,p1=%lf,p2=%lf,x0,x1,p0,p1); printf(na,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(f1f2) 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; return t; void gtd() double xN,gN,pN,t=0,f0,mod1=0,mod2=0,nanda=0; int i,k,n;printf(请输入函数的元数值n=); scanf(%d,&n); printf(n请输入初始值:n);for(i=0;ieps) p0=-g0;p1=-g1; k=0; while(1) t=hjfg(x,p);/*调用黄金分割法求t的值*/ printf(np1=%lf,p2=%lf,t=%lf,p0,p1,t); x0=x0+t*p0;x1=x1+t*p1; g0=2*x0; g1=50*x1; /*printf(nx1=%lf,x2=%lf,g1=%lf,g2=%lf,x0,x1,g0,g1);*/ mod2=sqrt(pow(g0,2)+pow(g1,2); /*求梯度的长度*/ if(mod2=eps) break; else if(k+1=n) g0=2*x0; g1=50*x1; p0=-g0; p1=-g1; k=0; else nanda=pow(mod2,2)/pow(mod1,2); printf(nnanda=%lf,mod=%lf,nanda,mod2); p0=-g0+nanda*p0; p1=-g1+nanda*p1; mod1=mod2; k+; printf(n-); printf(n最优解为x1=%lf,x2=%lf,x0,x1); printf(n最终的函数值为%lf,f(x,g,t); main() gtd(); 附录10#include#include#includeint const n=2;/正定二次函数的自变量个数double fun(double xn,double f_xsn+n+1+(n-1)*n/2);/输入变量为函数自变量初值void Q_fun(double f_xsn+n+1+(n-1)*n/2,double Qnn+1);void D_fun(double xn,double Qnn+1,double g0n);int H(double g0n,double c);/判别准则：返回1结束，返回0继续迭代void abc(double xn,double pn,double f_xsn+n+1+(n-1)*n/2,double t3); double t_c(double t3);/二次函数一阶导为零计算t的值，t0void main()double f_xsn+n+1+(n-1)*n/2=4,1,-40,-12,136,0;double xn=8