最优化方法大作业
最优化方法大作业答案
1.用薄钢板制造一体积5m 3，长度不小于4m ，无上盖的货箱，要求钢板耗量最小。
确定货箱的长x 1、宽x 2和高x 3。
试列出问题的数学模型。
解：min 32312122x x x x x x z ++= s.t 5321=x x x 41≥x 0,,321≥x x x2.将下面的线性规划问题表示为标准型并用单纯形法求解max f=x 1+2x 2+x 3s ．t ．2x 1+x 2-x 3≤2 -2x 1+x 2-5x 3≥-6 4x 1+x 2+x 3≤6 x i ≥0 i=1，2，3 解：先化标准形：Min 321x x x z -+=224321=+-+x x x x 6525321=++-x x x x646321=+++x x x x列成表格：121610011460105122001112-----可见此表已具备1°，2°，3°三个特点，可采用单纯形法。
首先从底行中选元素-1，由2/2,6/2,6/4最小者决定选第一行第一列的元素2，标以记号，迭代一次得121210231040116201002121211--------再从底行中选元素-2/3，和第二列正元素1/2，迭代一次得12123230210231040116201002121211-------再从底行中选元素-3，和第二列正元素2，迭代一次得4233410120280114042001112---再迭代一次得1023021062210231010213000421021013--选取最优解：01=x 42=x 23=x3. 试用DFP 变尺度法求解下列无约束优化问题。
min f （X ）=4（x 1-5）2+（x 2-6）2取初始点X=（8，9）T ，梯度精度ε=0.01。
解：取IH=0,初始点()TX 9,8=2221)6()5(4)(-+-=x x x f⎥⎦⎤⎢⎣⎡--=∇122408)(21x x x f⎪⎪⎭⎫⎝⎛=∇624)()0(xfTx f d )6,24()()0()0(--=-∇=)0(0)0()1(dxxα+=T)69,248(00αα--=])669()5248(4min[)(min 2020)0(0)0(--+--⨯=+αααdxf 0)6()63(2)24()2458(8)(00)0(0)0(=-⨯-+-⨯--=+ααααd d xdf13077.0130170≈=α⎪⎪⎭⎫⎝⎛=⎪⎪⎭⎫ ⎝⎛--⨯+⎪⎪⎭⎫ ⎝⎛=21538.886153.462413077.098)1(x⎪⎪⎭⎫⎝⎛-=∇43077.410784.1)()1(xf进行第二次迭代：⎥⎦⎤⎢⎣⎡--=-=78463.013848.31)0()1(xxδ⎥⎦⎤⎢⎣⎡--=∇-∇=56924.110783.25)()(1)0()1(xf xf γ101011011101γγγγγδδδH HH H H TTTT-+=03172.8011=γδT86614.6321101==γγγγH T⎥⎦⎤⎢⎣⎡=61561.046249.246249.285005.911Tδδ⎥⎦⎤⎢⎣⎡==46249.240022.3940022.3940363.630110110TTHH γγγγ所以：⎪⎪⎭⎫⎝⎛--=0038.103149.003149.012695.01H⎪⎪⎭⎫⎝⎛-⨯⎪⎪⎭⎫⎝⎛---=∇-=43076.410784.10038.103149.003149.012695.0)()1(1)1(xf H d⎪⎪⎭⎫⎝⎛-=48248.428018.0令 )1(1)1()2(dx x α+=利用)()1()1(=+ααd dxdf ，求得49423.01=α，所以⎪⎪⎭⎫⎝⎛-+⎪⎪⎭⎫⎝⎛=+=21538.213848.021538.886152.449423.0)1()1()2(dxx⎪⎪⎭⎫ ⎝⎛=65因)()2(=∇xf ，于是停，)2(x 即为最优解。
最优化方法大作业-算法源程序-0.618法、抛物线法、共轭梯度法
最优化方法程序作业专业：油气储运工程班级：姓名：学号：一、开发工具该应用程序采用的开发工具是Visual studio 2005，编程语言使用的是C#。
二、程序代码（1）0.618法和抛物线法求ψ(t)=sint，初始点t0=1①主程序：using System;using System.Data;using System.Configuration;using System.Web;using System.Web.Security;using System.Web.UI;using System.Web.UI.WebControls;using System.Web.UI.WebControls.WebParts;using System.Web.UI.HtmlControls;using System.Collections;public partial class_Default : System.Web.UI.Page{JiSuan JS = new JiSuan();protected void Button1_Click(object sender, EventArgs e){double xx0 = 0.01;//步长double tt0 = 1, pp = 2, qq = 0.5;ArrayList list1 = new ArrayList();list1 = (ArrayList)JS.SuoJian(xx0, tt0, pp, qq);//调用倍增半减函数double aa = double.Parse(list1[0].ToString());double bb = double.Parse(list1[1].ToString());txtShangxian.Text = bb.ToString();//在页面上显示极小区间txtXiaxian.Text = aa.ToString();ArrayList list2 = new ArrayList();list2 = (ArrayList)JS.JiXiao1(aa, bb);//调用0.618法函数double jixiao1 = double.Parse(list2[0].ToString());double fjixiao1 = double.Parse(list2[1].ToString());txtJixiao1.Text = jixiao1.ToString();//在页面上显示极小点txtFjixiao1.Text = fjixiao1.ToString();//在页面上显示极小点处的函数值ArrayList list3 = new ArrayList();list3 = (ArrayList)JS.JiXiao2(aa, bb);//调用抛物线法函数double jixiao2 = double.Parse(list3[0].ToString());double fjixiao2 = double.Parse(list3[1].ToString());txtJixiao2.Text = jixiao2.ToString();//在页面上显示极小点txtFjixiao2.Text = fjixiao2.ToString();//在页面上显示极小点处的函数值 }}②各个子函数的程序：using System;using System.Data;using System.Configuration;using System.Web;using System.Web.Security;using System.Web.UI;using System.Web.UI.WebControls;using System.Web.UI.WebControls.WebParts;using System.Web.UI.HtmlControls;using System.Collections;///<summary>/// JiSuan 的摘要说明///</summary>public class JiSuan{public JiSuan(){//// TODO: 在此处添加构造函数逻辑//}//目标函数public double f(double z){double zz = Math.Sin(z);return zz;}//倍增半减函数public ArrayList SuoJian(double x0, double t0, double p, double q){double t1;double f0, f1;double x1, t2;double f2;double a = 0, b = 0, c = 0;double fa, fc, fb;double t3, f3;bool biaozhi1 = true;//设置是否循环的标志bool biaozhi2 = true;t1 = t0 + x0;f0 = f(t0);//调用目标函数f1 = f(t1);if (f0 > f1){while (biaozhi1){x1 = p * x0;t2 = t1 + x1;f2 = f(t2);//调用目标函数if (f1 > f2){t0 = t1;t1 = t2;f0 = f1;f1 = f2;continue;}else{a = t0;c = t1;b = t2;fa = f0;fc = f1;fb = f2;break;}}}else{t3 = t1;f3 = f1;t1 = t0;f1 = f0;t0 = t3;f0 = f3;while (biaozhi2){x1 = q * x0;t2 = t1 - x1;f2 = f(t2);//调用目标函数if (f1 > f2){t0 = t1;t1 = t2;f0 = f1;f1 = f2;continue;}else{a = t2;c = t1;b = t0;fa = f2;fc = f1;fb = f0;break;}}}ArrayList list = new ArrayList();list.Add(a);list.Add(b);return list;}//0.618法函数public ArrayList JiXiao1(double a, double b) {double jd = 0.0001;//精度double p = 0.618;double t1, t2;double f1, f2;double jixiao=0;//极小点bool biaozhi1 = true;//设置是否循环标志bool biaozhi2 = true;while (biaozhi1){t1 = a + (1 - p)*(b - a);t2 = a + p * (b - a);f1 = f(t1);//调用目标函数f2 = f(t2);while (biaozhi2){if (f1 < f2){b = t2;t2 = t1;f2 = f1;t1 = a + (1 - p)*(b - a); f1 = f(t1);//调用目标函数if (Math.Abs(b - a) > jd) {continue;}else{biaozhi1 = false;break;}}else if (f1 == f2){a = t1;b = t2;if (Math.Abs(b - a) > jd) {break;}else{biaozhi1 = false;break;}}else if (f1 > f2){a = t1;t1 = t2;f1 = f2;t2 = a + p * (b - a);f2 = f(t2);//调用目标函数if (Math.Abs(b - a) > jd){continue;}else{biaozhi1 = false;break;}}}}jixiao = (a + b) / 2;double fjixiao = f(jixiao);//调用目标函数ArrayList list = new ArrayList();list.Add(jixiao);list.Add(fjixiao);return list;}//抛物线法函数public ArrayList JiXiao2(double a, double b){double jd = 0.0001;//精度double c = a + (b - a) / 2;//将c的值设为a、b点的中点double fa, fb, fc, c1, c2;double jixiao = 0;//极小点double ta, fta;double fac;bool biaozhi1 = true;//设置是否循环标志while (biaozhi1){fa = f(a);//调用目标函数fb = f(b);fc = f(c);c1 = (fb - fa) / (b - a);c2 = ((fc - fa) / (c - a) - c1) / (c - b);if (c2 == 0){jixiao = c;break;}else{ta = 0.5 * (a + b - (c1 / c2));fta = f(ta);//调用目标函数if (Math.Abs(b - a) <= jd){jixiao = ta;break;}else{if (fc > fta){if (c > ta){b = c;fb = fc;c = ta;fc = fta;continue;}else{a = c;fa = fc;c = ta;fc = fta;continue;}}else if (fc == fta){if (c > ta){a = ta;b = c;c = (c + ta) / 2;fa = fta;fb = fc;fc = f(c);//调用目标函数continue;}else if (c == ta){fac = f((a + c) / 2);//调用目标函数if (fac < fc){c = (a + c) / 2;fc = f(c);//调用目标函数b = ta;fb = fta;continue;}else{a = (a + c) / 2;fa = f(a);//调用目标函数continue;}}else if (c < ta){a = c;c = (c + ta) / 2;b = ta;fa = fc;fb = fta;fc = f(c);//调用目标函数continue;}}else if (fc < fta){if (c > ta){a = ta;fa = fta;continue;}else{b = ta;fb = fta;continue;}}}}}double fjixiao = f(jixiao);//调用目标函数ArrayList list = new ArrayList();list.Add(jixiao);list.Add(fjixiao);return list;}}（2）共轭梯度法求函数极小点：f(x)=1.5x12+0.5x22-x1x2-2x1；初始点为（-2，4）①主程序：using System;using System.Data;using System.Configuration;using System.Web;using System.Web.Security;using System.Web.UI;using System.Web.UI.WebControls;using System.Web.UI.WebControls.WebParts;using System.Web.UI.HtmlControls;using System.Collections;public partial class_Default : System.Web.UI.Page{JiSuan JS = new JiSuan();protected void Button1_Click(object sender, EventArgs e){//共轭梯度法求极小值double[] x0 = new double[2];//定义二维数组double[] x = new double[2];x0[0] = -2;x0[1] = 4;double x00 = -2;double x01 = 4;double jd = 0.0001;//精度bool biaozhi1 = true;//设置是否循环的标志bool biaozhi2 = true;double y;//定义关于x的函数值double[] g = new double[2];//定义函数在点x处的梯度值double[] p = new double[2];double ggm, ggo=0;double jixiao;double x1=0, x2=0, fy=0;int i;for (int j = 0; j < 2; j++){x[j] = x0[j];}while (biaozhi1){i = 0;while (biaozhi2){y = JS.f(x[0], x[1]);//调用目标函数g[0] = JS.g0(x[0], x[1]);//调用梯度函数，求在点x处的梯度值 g[1] = JS.g1(x[0], x[1]);ggm = g[0] * g[0] + g[1] * g[1];if (ggm <= jd){x1 = x[0];x2 = x[1];fy = y;biaozhi1 = false;break;}else{if (i == 0){p[0] = -g[0];p[1] = -g[1];}else{p[0] = -g[0] + (ggm / ggo) * p[0];p[1] = -g[1] + (ggm / ggo) * p[1];}//调用0.618法子函数，找出极小点jixiao = JS.JiXiao1(x[0], x[1], p[0], p[1], x00, x01); x[0] = x[0] + jixiao * p[0];x[1] = x[1] + jixiao * p[1];ggo = ggm;i++;if (i >= 2){break;}else{continue;}}}}txtx1.Text = x1.ToString();txtx2.Text = x2.ToString();txty.Text = fy.ToString();}}②子函数程序：using System;using System.Data;using System.Configuration;using System.Web;using System.Web.Security;using System.Web.UI;using System.Web.UI.WebControls;using System.Web.UI.WebControls.WebParts;using System.Web.UI.HtmlControls;using System.Collections;///<summary>/// JiSuan 的摘要说明///</summary>public class JiSuan{public JiSuan(){//// TODO: 在此处添加构造函数逻辑//}//目标函数public double f(double z1,double z2){double zz = 1.5 * z1 * z1 + 0.5 * z2 * z2 - z1 * z2 - 2 * z1;return zz;}//求梯度值函数public double g0(double z1, double z2){double zz = 3 * z1 - z2 - 2;return zz;}public double g1(double z1, double z2){double zz = z2 - z1;return zz;}//倍增半减函数public ArrayList SuoJian(double x0, double t0, double p, double q, double xx0, double xx1, double p0, double p1,double x00,double x01){double t1;double f0, f1;double x1, t2;double f2;double a = 0, b = 0, c = 0;double fa, fc, fb;double t3, f3;bool biaozhi1 = true;//设置是否循环的标志bool biaozhi2 = true;double x_0t1, x_1t1, x_0t2, x_1t2;x_0t1 = xx0 + t1 * p0;x_1t1 = xx1 + t1 * p1;f0 = f(x00, x01);//调用目标函数f1 = f(x_0t1, x_1t1);if (f0 > f1){while (biaozhi1){x1 = p * x0;t2 = t1 + x1;x_0t2 = x_0t1 + t2 * p0;x_1t2 = x_1t1 + t2 * p1;f2 = f(x_0t2, x_1t2);//调用目标函数if (f1 > f2){t0 = t1;t1 = t2;f0 = f1;f1 = f2;continue;}else{a = t0;c = t1;b = t2;fa = f0;fc = f1;fb = f2;break;}}}else{t3 = t1;f3 = f1;t1 = t0;f1 = f0;t0 = t3;while (biaozhi2){x1 = q * x0;t2 = t1 - x1;x_0t2 = x_0t1 + t2 * p0;x_1t2 = x_1t1 + t2 * p1;f2 = f(x_0t2, x_1t2);//调用目标函数if (f1 > f2){t0 = t1;t1 = t2;f0 = f1;f1 = f2;continue;}else{a = t2;c = t1;b = t0;fa = f2;fc = f1;fb = f0;break;}}}ArrayList list = new ArrayList();list.Add(a);list.Add(b);return list;}//0.618法函数public double JiXiao1(double x0, double x1, double p0, double p1,double x00,double x01){double jd = 0.0001;//精度double p = 0.618;double t1, t2;double jixiao = 0;//极小点bool biaozhi1 = true;//设置是否循环标志bool biaozhi2 = true;double xx0 = 0.01;//步长double tt0 = 0, pp = 2, qq = 0.5;double x_0t1, x_1t1, x_0t2, x_1t2;ArrayList list1 = new ArrayList();//调用倍增半减函数获取极小区间list1 = (ArrayList)SuoJian(xx0, tt0, pp, qq, x0, x1, p0, p1, x00, x01);double a = double.Parse(list1[0].ToString());double b = double.Parse(list1[1].ToString());while (biaozhi1){t1 = a + (1 - p) * (b - a);t2 = a + p * (b - a);x_0t1 = x00 + t1 * p0;x_1t1 = x01 + t1 * p1;x_0t2 = x_0t1 + t2 * p0;x_1t2 = x_1t1 + t2 * p1;f1 = f(x_0t1, x_1t1);//调用目标函数f2 = f(x_0t2, x_1t2);while (biaozhi2){if (f1 < f2){b = t2;t2 = t1;f2 = f1;t1 = a + (1 - p) * (b - a);x_0t1 = x00 + t1 * p0;x_1t1 = x01 + t1 * p1;f1 = f(x_0t1, x_1t1);//调用目标函数if (Math.Abs(b - a) > jd){continue;}else{biaozhi1 = false;break;}}else if (f1 == f2){a = t1;b = t2;if (Math.Abs(b - a) > jd){break;}else{biaozhi1 = false;break;}}else if (f1 > f2){a = t1;t1 = t2;f1 = f2;t2 = a + p * (b - a);x_0t2 = x_0t1 + t2 * p0;x_1t2 = x_1t1 + t2 * p1;f2 = f(x_0t2, x_1t2);//调用目标函数if (Math.Abs(b - a) > jd){continue;}else{biaozhi1 = false;break;}}}}jixiao = (a + b) / 2;return jixiao;}}三、程序运行界面及结果（1）0.618法和抛物线法求ψ(t)=sint，初始点t0=1（2）共轭梯度法求函数极小点：f(x)=1.5x12+0.5x22-x1x2-2x1；初始点为（-2，4）。
优化方法大作业1
优化方法大作业一、Wolfe-Powell法利用MATLAB软件编写，其中初始值x0=-5；其他参数按照已知条件来取。
当b分别等于8、9、10时，均得到如下结果：而当初值x0变化时，则结果变化比较大，如将x0取-6，则计算结果如下：通过比较可以看出，b的取值对计算结果的影响较小，而初始值x0的取值则对结果影响很大。
从中也表明Wolfe-Powell准则的收敛条件比较弱，容易出现当函数还没取极小值而迭代循环已结束的情况。
具体代码见附录1.二、无约束优化1.DFP法精度ε = 0.02，确定λ使用的是一维非精确搜索算法(直接法，Wolfe-Powell准则)。
结果如下图所示：2. BFGS方法同上一种方法，精度ε = 0. 02，确定λ使用的是一维非精确搜索算法(直接法，Wolfe-Powell准则)。
结果如下图所示：比较两种方法，查看每一步的函数值，得出如下结果。
图1：DFP与BFGS法结果对比图通过图1与表1的比较，BFGS比DFP法最初收敛得更快，但是DFP法比BFGS法的最终结果更好。
DFP与BFGS法的代码分别见附录2、3。
三、Rockafellar乘子法文件myfun.m：function y=myfun(x)y=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3);文件mycon.m：function [c,ceq]=mycon(x)c(1)=(-x(1));c(2)=(-x(2));c(3)=(-x(3));ceq(1)=x(1)^2+x(2)^2+x(3)^2-25;ceq(2)=8*x(1)+14*x(2)+7*x(3)-56文件Q3.m：clearclcx0=[2 2 2]';[x,fval,exitflag,output]=fmincon(@myfun,x0,[],[],[],[],[],[],@mycon);附录1：%% Wolfe-Powell 准则方法clear;tic;c1=0.1;c2=0.65;a=0;b=8;syms x;f=x*x-2*x+7; %函数G0=jacobian(f,x);%求梯度x0=-6; %初始取x=-5% 由于只有一个未知数，默认S=1lambda = 1;x=x0;for k=1:1:1000f0 = eval(f); %求出函数值grad1= eval(G0); %求出梯度值figure1(k) = f0; % 用于绘出迭代过程中函数值变化x= x0+lambda;f1 = eval(f); %求出函数值grad2 = eval(G0); %求出梯度值value1 = f0 - f1 + c1* lambda * grad1 ;value2 = grad2 -c2*grad1 ;if value1 >=-1e-6 && value2 >=-1e-6disp('The variable matrix is : ');disp(x);fun=eval(f);fprintf('The minimum value of function is %f \n',fun);break; %如果满足两个条件，则退出循环。
实用最优化方法编程大作业
实用最优化方法编程大作业班级：姓名：指导老师：学号：大连理工大学2015年11月27日版本号：MATLAB 7.11.0 (R2010b)【文件名WP.m】function x = WP(f,x0,var,s0,eps)clcsyms x1 x2 ;c1=0.1;c2=0.5;b=inf;lambda=1;ifnargin==4eps=1.0e-6;endgradf=jacobian(f,var);g0=subs(gradf,var,x0);f0=subs(f,var,x0);gs=g0*s0';a=0;j=0;while j<1000new_x=x0+lambda*s0;new_f=subs(f,var,new_x);left=f0-new_f;new_g=subs(gradf,var,new_x);new_gs=new_g*s0';right=-1*c1*lambda*gs;if left<right %不满足第一个条件j=j+1;b=lambda;lambda=0.5*(lambda+a);else %满足第一个条件left2=new_gs;right2=c2*g0;if left2<right2 %不满足第二个条件j=j+1;a=lambda;lambda=min(2*lambda,0.5*(lambda+b));elsex=lambda;break;endendend在Command Window 输入：symsx1 x2WP(100*(x2-x1^2)^2+(1-x1)^2,[-1,1],[x1 x2],[1 1])程序运行后可得出结果：ans=0.0039【文件名minGETD.m】function [x,minf] = minGETD(f,x0,var,eps)clcsyms x1 x2 x3format long;n=length(x0);ifnargin==3eps=1.0e-3;endx0=x0';syms lambda;gradf=jacobian(f,var);g=subs(gradf,var,x0);s=-g';k=0;while 1tol=norm(double(g));iftol<=epsx=x0;break;endx1=x0+lambda*s;f1=subs(f,var,x1);dy1=diff(f1);lambda0=solve(dy1);x1=x0+lambda0*s;g1=subs(gradf,var,x1)tol=norm(double(g1))iftol<=epsx=x1;break;endif k+1==nx0=x1;continue;elsemiu=dot(g1,g1)/dot(g,g)s=-g1'+miu*s;k=k+1;x0=x1;endendx在Command Window 输入：symsx1 x2 x3x0=[0 0 0];var=[x1 x2 x3];f=x1^2-2*x1*x2+2*x2^2+x3^2-x2*x3+2*x1+3*x2-x3;minGETD(f,x0,var,eps)程序运行后可得出结果：x=[-236894/59711,-178563/59711,-59465/59711]可认为最终解为[-4,-3,-1]。
最优化方法与应用大作业(一)最速下降法
最优化方法与应用大作业（一）
---最速下降法部分：
1.问题描述:
用梯度下降法求解以下优化问题
min f(x)=(x1+10*x2)^2+5(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4
2.编程感想：
该算法需要计算Hesse矩阵，C语言在向量运算时没有Matlab方便，所以手工完成了理论计算，再输入，破坏了程序的移植性。
同时，实验表明当初始值离理想点较远且精度要求较高时，最速下降法的收敛速率极慢，迭代几乎不可能完成，这对初值的选取提出了一定限制。
3.结果分析：
编译界面（Mac os X，Xcode环境）
输入参数（设定为(0.1,0.2,0.3,0.4)）：
结果（此处列出每次迭代结果）。
明显的看到，最速下降法的收敛较慢，最终结果接近理论值（F（0，0，0，0）=0）所以该结果可以满意。
4.算法代码见下页
西安电子科技大学电子工程学院020951
李骏昊02095005。
最优化方法大作业模板
命题人：审核人：大作业学期：至学年度第学期课程：最优化方法课程代号：签到序号：使用班级：姓名：学号：题号一二三四五六七八九十总分得分一、（目标1）请从以下6种算法中任选一种，说明算法的来源、定义、基本思想和优缺点，并给出算法步骤（包含算法流程图）和例子（包含程序与运算结果）。
①禁忌搜索算法；②模拟退火算法；③遗传算法；④神经网络算法；⑤粒子群算法；⑥蚁群算法。
二、（目标1）某工厂生产甲、乙两种产品，已知生产这两种产品需要消耗三种材料A 、B 和C ，其中生产过程中材料的单位产品消耗量和总量，以及单位产品的利润如下表所示。
该如何配置安排生产计划，使得工厂所获得的利润最大？材料甲乙资源总量材料A （Kg ）3265材料B （Kg ）2140材料C （Kg ）0375单位利润（元/件）15002500-（1）要保证工厂利润的最大化，写出相应的生产计划数学模型；（2）根据对偶理论，直接写出该线性规划的对偶问题；（3）采用单纯形表法对该该线性规划问题进行求解，写出详细的计算过程；（4）采用Matlab 软件对该线性规划问题进行求解，写出完整的源程序，并给出程序运行结果；（5）讨论当材料B 的资源总量发生变化时，该线性规划问题的最优解会如何变化？课程目标目标1……题号一、二、三、四、五……分值20、25、20、20、15……得分得分三、（目标1）求解下列无约束非线性规划问题（1）采用黄金分割法求解：min 4()24f x x x =++。
初始区间为[-1.0]，精度为ε=10-4。
（要求：采用黄金分割法进行Matlab 编程求解，写出源程序，并给出运行结果，列出迭代过程的数据表格）（2）采用阻尼牛顿法求解：222121212min (,)4f x x x x x x =+-。
分别取两个初始点：x A =(1,1)T ，x B =(3,4)T 。
（要求：采用阻尼牛顿法进行Matlab 编程求解，并给出运行结果，列出迭代过程的数据表格）四、（目标1）求解下列约束非线性规划问题：22112212121212min ()23532..00f x x x x x x x x x x x s t x x =-+--+≤⎧⎪-≤⎪⎨≥⎪⎪≥⎩（1）采用罚函数法进行求解，需写出具体计算过程；（2）采用二次规划方法进行求解，需写出具体计算过程，并进行MATLAB 编程，写出源程序和运算结果；五、（目标1）（1）某商店在未来的4个月里，准备利用它的一个仓库来专门经营某种商品，仓库的最大容量为1000单位，而且该商店每月只能出卖仓库现有的货。
最优化方法课程设计-最优化大作业-用优化算法求解函数最值问题
最优化方法大作业---------用优化算法求解函数最值问题摘要最优化(optimization) 是应用数学的重要研究领域.它是研究在给定约束之下如何寻求某些因素(的量)，以使某一(或某些)指标达到最优的一些学科的总称。
最优化问题一般包括最小化问题和最大化问题，而最大化问题可以通过简单的转化使之成最最小化问题。
最小化问题分为两类，即约束最小化和无约束最小化问题。
在此报告中，前两个问题属于无约束最小化问题的求解，报告中分别使用了“牛顿法”和“共轭梯度法”。
后两个问题属于有约束最小化问题的求解，报告中分别用“外点法”和“内点法”求解。
虽然命名不一样，其实质都是构造“惩罚函数”或者“障碍函数”，通过拉格朗日乘子法将有约束问题转化为无约束问题进行求解。
再此报告中，“外点法”和“内点法”分别用了直接求导和调用“牛顿法”来求解无约束优化问题。
在此实验中，用“共轭梯度法”对“牛顿法”所解函数进行求解时出现错误，报告中另取一函数用“共轭梯度法”求解得到正确的结果。
此实验中所有的函数其理论值都是显见的，分析计算结果可知程序正确，所求结果误差处于可接受范围内。
报告中对所用到的四种方法在其使用以前都有理论说明，对“外点法”中惩罚函数和“内点法”中障碍函数的选择也有相应的说明，另外，对此次试验中的收获也在报告的三部分给出。
本报告中所用程序代码一律用MATLAB编写。
【关键字】函数最优化牛顿法共轭梯度法内点法外点法 MATLAB一，问题描述1，分别用共轭梯度发法和牛顿法来求解一下优化问题()()()()()441432243221102510min x x x x x x x x x f -+-+-++=2， 分别用外点法和内点发求解一下优化问题⎩⎨⎧≥-++01.min 212231x x t s x x二、问题求解1.1 用牛顿法求解()()()()()441432243221102510min x x x x x x x x x f -+-+-++=1.1.1问题分析：取步长为1而沿着牛顿方向迭代的方法称为牛顿法，在牛顿法中，初始点的取值随意，在以后的每次迭代中，()[]()k k k k x f x f x x ∇∇-=-+121，直到终止条件成立时停止。
最优化方法大作业
单位代码03学号《最优化方法》课程实践完成时间：2015年5月30日星期六选择题目：题目一使用优化软件，编写重要算法的程序1.第一大题：（1）学习最优流量工程问题，nonsmooth_MCFP.pdf（2）问题重述：Figure 1一个简单的网络拓扑和流量需求如Figure 1所示，网络有7 个节点，13 条弧，每条弧的容量是5 个单位. 此外有四个需求量均为4个单位的源－目的对（M=4)，具体的源节点、目的节点信息如图所示. 这里为了简单，省去了未用到的弧，此外弧上的数字表示弧的编号。
（3）极小化MAU设定变量x，为531⨯的向量，其中(53)x即为变量z。
使用linprog 函数求解极小化问题得到x。
之前确定三个约束条件。
1、Ax b≤，其中A为1353⨯的矩阵，b为131⨯的向量。
2、eq eq x b A =，其中eq A 为2853⨯的矩阵，eq b 为281⨯的向量。
3、x lb ≥，其中lb 为153⨯的向量 编程计算后得到结果如下：（4） 极小化FT 成本函数设定变量x ，为651⨯的向量，其中(53:65)x 即为变量l z 。
使用linprog 函数求解极小化问题得到x 。
之前确定三个约束条件。
1、Ax b ≤，其中A 为7865⨯的矩阵，b 为781⨯的向量。
2、eq eq x b A =，其中eq A 为2865⨯的矩阵，eq b 为281⨯的向量。
3、x lb ≥，其中lb 为165⨯的向量 编程计算后得到结果如下：2. 第二大题： 2.1. 习题5.6 2.1.1. 问题分析问题2112212()(101810)/241513q x x x x x x x =-++-+ 通过matlab 画出其等高线为：2.1.2. 最速下降法最速下降法中，取值：k k p g =-==()()k k k kk T k k T k g p g g p Gp g Ggα- x1x 2等高线-224681012(1)()k x k x k p α+=+2.1.3. 算法流程图如下图所示：2.1.4. 初始值（0,0）编程运行结构为：收敛过程曲线为：2.1.5. 初始值（-0.4,0）编程运行结构为：收敛过程曲线为：x1x 2等高线-2246810122.1.6. 初始值（10,0）编程运行结构为：收敛过程曲线为：x1x 2-2246810122.1.7. 初始值（11,0）编程运行结构为：收敛过程曲线为：x1x 2-2246810122.2. 习题5.7 2.2.1. 问题分析问题()94ln(7)f x x x =--497g x =-- 24(7)G x =- Matlab 画出在区间(7 10)的函数、一阶导数、二阶导数的变化曲线为x1x 2-22468101277.588.599.510707274767880828486xf函数变化曲线77.588.599.510-35-30-25-20-15-10-50510xg一阶导数g 变化曲线2.2.2. 牛顿法牛顿法中，取值：k k k G s g =-1k k k s xx +=+其中，如果G 不是半正定，则采用修正牛顿法(+)k k k G I s g λ=-77.588.599.51050100150200250300350400xg二阶导数G 变化曲线2.2.3.算法流程图如下图所示：2.2.4.初始值7.40编程运行结构为：收敛过程曲线为：2.2.5. 初始值7.20编程运行结构为：收敛过程曲线为：7.397.47.417.427.437.447.457.4670.24570.2570.25570.2670.265xf2.2.6. 初始值7.01编程运行结构为：收敛过程曲线为：7.17.27.37.47.57.67.770.270.470.670.87171.271.471.6xf2.2.7. 初始值7.80编程运行结构为：收敛过程曲线为：6.856.9 6.9577.057.17.157.27.257.37.35727476788082xf2.2.8. 初始值7.88编程运行结构为：收敛过程曲线为：7.17.27.37.47.57.67.77.87.97070.57171.572xf2.2.9. 分析函数在区间（7，7.8888）内是凸函数，G 恒大于零，所以单纯牛顿法保证收敛。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
发动机空燃比控制器引言：我主要从事自动化相关研究。
这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。
发动机空燃比控制器设计中的最优化问题AFR =afm m && (1)空燃比由方程（1）定义，在发动机运行过程中如果控制AFR 稳定在14.7可以获得最好的动力性能和排放性能。
如果假设进入气缸的空气流量am &可以由相关单元检测得到，则可以通过控制进入气缸的燃油流量f m &来实现空燃比的精确控制。
由于实际发动机的燃油喷嘴并不是直接对气缸喷燃油，而是通过进气歧管喷燃油，这么做会在进气歧管壁上液化形成油膜，因此不仅是喷嘴喷出的未液化部分燃油会进入气缸，油膜蒸发部分燃油也会进入气缸，如方程（2）。
这样如何更好的喷射燃油成为了一个问题。
1110101122211ττττ⎡⎤⎡⎤-⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎣⎦⎢⎥⎣⎦f ff v X x x u x x X x y =x && (2)其中12、,==ff fv x m x m &&=f y m &,=fi u m &这里面，表示油膜蒸发量ff m &、fvm &表示为液化部分燃油、fim &表示喷嘴喷射的燃油，在τf 、τv 、X 都已知的情况下，由现代控制理论知识，根据系统的增广状态空间模型方程（3）00000011011011114.70ττττ⎡⎤⎡⎤-⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥=-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤=⎢⎥⎣⎦⎣⎦ff v v a X X u +q q m y q x x x &&& (3)其中()014.7⎰taq =y -m&。
由极点配置方法，只要设计控制器方程（4），就可以使得y 无差的跟踪阶跃输入，那么y 也能较好的跟踪AFR *am /&。
12-- u =K q K x (4)这里面的12、K K 确定，可由主导极点概念降维成两个参数12C ,C ，虽然都是最终稳态无差，但是目标是使得瞬态过程中y 和阶跃输入y r 的差异尽可能的小。
所以原问题转化成了一个关于参数的12C ,C 的，以()()2=-⎰W r c y y 为目标函数的优化问题。
这里面由于无法写出参数与目标函数的显示表达式，从而也就无法使用和导数有关的优化方法。
所以我选择用坐标轮换法，由于坐标轮换法，局部寻优性能好，且不需要事先知道函数的导数。
在matlab7.10.0上对其进行仿真试验（程序见附录）。
假设进入气缸的空气质量流量已知，并且已经适当量化，假设状态值已经由状态观测器测量得到，在τf = 1s 、τv = 60ms 、X = 0.5时按上述方法设计控制器。
坐标轮换法的精度10.051ε⎡⎤=⎣⎦,a ,b其中，11001010，''⎡⎤⎡⎤==⎣⎦⎣⎦a ,b ,为12'⎡⎤⎣⎦c ,c 的初次搜索区间；基于阶跃响应的系统性能指标()()2=-⎰W r c y y ，通过坐标轮换法得到最优参数1 4.0542,0.9544,2**=-=c c 得到最优极点1 4.0542+0.9544, 4.05420.92λλ**=-=--j 544j ，将最优极点带入爱克曼公式得到反馈增益矩阵9.7820,4.4379,43.3499⎡⎤-⎣⎦=K ，进一步得到控制器：()9.7820,4.437943.349914.7⎡⎤⎡⎤⎢⎥=---*-⎣⎦⎢⎥⎣⎦fffi a f fv m m mm m &&&&& （5） 画出最后的仿真跟踪效果图，如图1，可见坐标轮换法很好的完成了任务。
图1 fm &对14.7a m &的跟踪效果从而使得发动机空燃比的变化如图2。
图2 空燃比变化time(sec)time(sec)空燃比主函数tf = 1;tv = 0.06; X=0.5;k0 = 10000000;k = 0;A = [-1/tf 0;0 -1/tv];tf0 = 1;tv0 = 0.06;X0 = 0.5;A1 = [-1/tf0 0;0 -1/tv0];B1 = [X0/tf0 (1-X0)/tv0]';C = [1 1];D = 0;AA = [A1 zeros(2,1);C 0];BB = [B1 ;0];P = zblh();PJ = [-10*sqrt(P(1)*P(1)+P(2)*P(2)) (-P(1)+P(2)*1i) (-P(1)-P(2)*1i)];K = acker(AA, BB, J)K1 = [K(1) K(2)];K2 = K(3);-10*sqrt(P(1)*P(1)+P(2)*P(2))Ac = [(A-B*K1) -B*K2;C 0];Bc = [0;0;-1];Cc = [C 0];Dc = 0;hold on;sys = ss(Ac,Bc,Cc,Dc);a1 = 0:0.01:3;b1 = 1+0*a1;a2 = 3.01:0.01:4;b2 = 1+(a2-3);a3 = 4.01:0.01:6;b3 = 2+0*a3;a4 = 6.01:0.01:10;b4 = 2+(a4-6)/2;a5 = 10.01:0.01:14;b5 = 4+0*a5;a6 = 14.01:0.01:17;b6 = 4-(a6-14);a7 = 17.01:0.01:20;b7 = 1+0*a7;t0 = [a1 a2 a3 a4 a5 a6 a7];u0 = [b1 b2 b3 b4 b5 b6 b7];x0 = [0.5 0.5 -0.1769]';[y0,t0,x] = lsim(sys,u0,t0,x0);x1 = [1 0 0 ]*x';x2 = [0 1 0]*x';x3 = [0 0 1]*x';x = [x1;x2];u = -K1*x-K2*x3;l = u0./y0';hold on;plot(t0,u0,'-');plot(t0,y0,':');xlabel('time(sec)')ylabel('Output')zblu（）坐标轮换法函数function [ R ] = zblh( )m = 0;n = 10;x2 = m + 0.618*(n-m);f2 = wch(x2,0.5);x1 = m + 0.382*(n-m);f1 = wch(x1,0.5);M = 0;N = 10;y2 = M + 0.618*(N-M);F2 = wch(1,y2);y1 = M + 0.382*(N-M);F1 = wch(1,y1);flag = 1;a = 0.5;b = 0.5;k1 = 0;k2 = 0;while(abs(abs(m-n)+abs(M-N))>0.05) if flag ~= 0k1 = k1 + 1;if f1 < f2n = x2;x2 = x1;f2 = f1;x1 = m + 0.382*(n-m);f1 = wch(x1,b);else if f1 == f2;m = x1;n = x2;x2 = m + 0.618*(n-m); f2 = wch(x2,b);x1 = m + 0.382*(n-m); f1 = wch(x1,b);elsem = x1;x1 = x2;f1 = f2;x2 = m + 0.618*(n-m); f2 = wch(x2,b);endenda = (m+n)/2;flag = 0;elsek2 = k2 +1;if F1 < F2N = y2;y2 = y1;F2 = F1;y1 = M + 0.382*(N-M);F1 = wch(a,y1);else if F1 == F2;M = y1;N = y2;y2 = M + 0.618*(N-M); F2 = wch(a,y2);y1 = M + 0.382*(N-M);F1 = wch(a,y1);elseM = y1;y1 = y2;F1 = F2;y2 = M + 0.618*(N-M);F2 = wch(a,y2);endendb = (M+N)/2;flag = 1;endendk1k2R = [a b]end目标函数function [ f ] = wch(a,b)tf = 1;tv = 0.06;X=0.5;k0 = 10000;k = 0;A = [-1/tf 0;0 -1/tv];tf0 = 1;A1 = [-1/tf0 0;0 -1/tv];B = [X/tf (1-X)/tv]';B1 = [X/tf0 (1-X)/tv]';C = [1 1];D = 0;AA = [A1 zeros(2,1);C 0];BB = [B1 ;0];J = [-10*sqrt(a*a+b*b) (-a+b*1i) (-a-b*1i)]; K = acker(AA, BB, J);K1 = [K(1) K(2)];K2 = K(3);Ac = [(A-B*K1) -B*K2;C 0];Bc = [0;0;-1];Cc = [C 0];Dc = 0;t0 = 0:0.01:4;yn = heaviside(t0);sys = ss(Ac,Bc,Cc,Dc);u0 = heaviside(t0);y0 = lsim(sys,u0,t0);for t0 = 1:400;kn = (yn(t0) - y0(t0))*(yn(t0) - y0(t0)); k = k+kn;endf = k;end。