《MATLAB解方程与最优化问题求解.ppt》由会员分享，可在线阅读，更多相关《MATLAB解方程与最优化问题求解.ppt（86页珍藏版）》请在三一办公上搜索。
1、MATLAB程序设计教程（第二版）,刘卫国 主编 中国水利水电出版社,第6章 MATLAB解方程与最优化问题求解,MATLAB线性方程组求解 MATLAB非线性方程数值求解 MATLAB常微分方程初值问题的数值解法 MATLAB最优化问题求解,背景经济与工程中的许多问题最后都可以转化为求解线性方程组许多数值计算问题（如样条函数、常微分方程数值解、差分方程等）的研究也往往归结为此类问题线性方程组的数值法直接法经过有限次算法运算求出精确解，最常用的是：高斯消元法、矩阵LU分解迭代法从初值出发，用递推的方法，给出近似解序列。常用的方法：雅可比迭代法、高斯-赛德尔迭代法。,列昂杰夫的“投入-产出”模型
2、：列昂惕夫用线性代数研究经济数学模型，1949年曾用计算机计算出了由美国统计局的25万条经济数据所组成的42个未知数42个方程的方程组。投入产出分析的方法基础包括：线性方程组和矩阵运算（静态模型）、微分方程和差分方程（动态模型）、电子计算机。,第三章 矩阵代数,设有n个经济部门，xi为部门i的总产出，cij为部门j单位产品对部门i产品的消耗，di为外部对部门i的需求，fj为部门j新创造的价值。分配平衡方程组消耗平衡方程组 i=1,2,n,第三章 矩阵代数,投入产出分析,令 C=（cij），X=(x1,xn),D=(d1,dn)，F=(f1,fn),则 X=CX+D令 A=EC，E为单位矩阵，则
3、 AX=DC称为直接消耗矩阵A称为列昂杰夫（Leontief）矩阵。,2023/11/8,7,第三章 矩阵代数,Y=1,1,1 B,Y表示各部门的总投入，称为投入向量。,新创造价值向量 F=X Y,B=C,B表示各部门间的投入产出关系，称为投入产出矩阵。,其中,常数项,所谓一般线性方程组是指具有形式：,由m个方程n个未知量的线性方程构成的方程组,矩阵形式：,6.1 线性方程组求解6.1.1 直接解法1利用左除运算符的直接解法对于线性方程组可以利用左除运算符“”求解：x=Ab；(1)若A为nn方阵，自动利用高斯消元法求解，若b是n1的列向量，则解x为n1的列向量,若b是nm的矩阵，可得到m个以A
4、为系数矩阵的线性方程组的数值解x(nm的矩阵),即x(:,j)=Ab(:,j),j=1,2,m.(2)当A不是方阵时，Ax=b称为欠定方程组或超定方程组，MATLAB会在最小二乘意义下解,例6-1 用直接解法求解下列线性方程组。,命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;%右边常数项应该是列向量x=Ab,2利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。优点：运
5、算速度快，节省存储空间,(1)LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。,MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：（1）L,U=lu(X)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。此时矩阵L往往不是下三角阵，但通过行交换后使之成为一个下三角阵（2）L,U,P=lu(X)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组A
6、x=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。,矩阵LU分解，设,命令如下：A=1,-1,1;5,-4,3;2,1,1;L,U=lu(A)%第（1）种格式LU=L*U%检验分解是否正确L,U,P=lu(A)LU=L*U%检验分解是否正确LU=inv(P)*L*U%检验分解是否正确,例6-2 用LU分解求解例6-1中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;L,U=lu(A);x=U(Lb)%或采用LU分解的第2种格式，命令如下：L,U,P=lu(A);x=U(LP*b),(2)QR分解对
7、矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：(1)Q,R=qr(X)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。(2)Q,R,E=qr(X)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb)。,矩阵QR分解，设,命令如下：A=1,-1,1;5,-4,3;2,7,10;Q,R=qr(A)%第（1）种格式QR=Q*R%检验分解是否正确Q,R,E=qr(A)QR
8、=Q*R%检验分解是否正确QR=Q*R/E%检验分解是否正确,例6-3 用QR分解求解例6-1中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;Q,R=qr(A);x=R(Qb)%或采用QR分解的第2种格式，命令如下：Q,R,E=qr(A);x=E*(R(Qb),(3)Cholesky分解如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即X=RR。,MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为：（
9、1）R=chol(X)：产生一个上三角阵R，使RR=X。若X为非对称正定，则输出一个出错信息。（2）R,p=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足RR=X(1:q,1:q)。实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R(Rb)。,矩阵Cholesky分解，设,命令如下：A=2,1,1;1,2,1;1,-1,3;R=chol(A)%第（1）种格式R*R%检验分解是否正确R,p=chol(A),p=0,则A是正定矩阵，对一个非
10、正定矩阵进行Cholesky分解，将得出错误信息，因而该方法可用来判断矩阵是否是正定矩阵,例6-4 用Cholesky分解求解例6-1中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;R=chol(A)?Error using=cholMatrix must be positive definite命令执行时，出现错误信息，说明A为非正定矩阵。,6.1.2 迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。,
11、1Jacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，即aii0(i=1,2,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是Jacobi迭代公式。如果序列x(k+1)收敛于x，则x必是方程Ax=b的解。,Jacobi迭代法的MATLAB函数文件Jacobi.m如下：function y,n=jacobi(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=e
12、ps x0=y;y=B*x0+f;n=n+1;end,例6-5 用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。,在命令中调用函数文件Jacobi.m，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6),2Gauss-Serdel迭代法在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为Gauss
13、-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。,Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function y,n=gauseidel(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y;y=G*x0+f;n=n+1;end,例6-6 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。,在命令中调用函数文件gauseidel.m，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=
14、9,7,6;x,n=gauseidel(A,b,0,0,0,1.0e-6),例6-7 分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。,命令如下：a=1,2,-2;1,1,1;2,2,1;b=9;7;6;x,n=jacobi(a,b,0;0;0)x,n=gauseidel(a,b,0;0;0),6.2 非线性方程数值求解6.2.1 单变量非线性方程求解 在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为：z=fzero(fname,x0,tol,trace)其中fname是待求根的函数文件名，x0为搜索的起点。一个
15、函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。,例6-8 求f(x)=x-10 x+2=0在x0=0.5附近的根。步骤如下：(1)建立函数文件funx.m。function fx=funx(x)fx=x-10.x+2;(2)调用fzero函数求根。z=fzero(funx,0.5)z=0.3758,6.2.2 非线性方程组的求解 对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：X=fsolve(fun
16、,X0,option)其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,例6-9 求下列非线性方程组在(0.5,0.5)附近的数值解。(1)建立函数文件myfun.m。f
17、unction q=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y);(2)在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。x=fsolve(myfun,0.5,0.5,optimset(Display,off)x=0.6354 0.3734,将求得的解代回原方程，可以检验结果是否正确，命令如下：q=myfun(x)q=1.0e-009*0.2375 0.2957 可见得到了较高精度的结果。,引例,2.讨论资金积累、国民收入、与人口增长的关系.（1）若国民平均
18、收入x与按人口平均资金积累y成正比，说明仅当总资金积累的相对增长率k大于人口的相对增长率r时，国民平均收入才是增长的.（2）作出k(x)和r(x)的示意图，分析人口激增会引起什么后果.,1.地中海鲨鱼问题：意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？,微分方程的解析解,结 果：u=tg(t+c),解 输入命令：dsolve(Du=1+u2,t),解 输
19、入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结 果 为:y=3e-2xsin（5x）,解 输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t);x=simple(x)%将x化简 y=simple(y)z=simple(z),结 果 为：x=c2e2t+c3e-t y=c1e-2t+c2e2t+c3e-t z=c1e-2t+c2e2t,结 果 为：x=C2*exp(2*t)+C3*exp(-t)y=C2*exp(2*t)+C3*exp(-t)+exp(-2*t)*C1
20、 z=C2*exp(2*t)+exp(-2*t)*C1,微分方程的数值解,常微分方程数值解的定义,在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。,因此，研究常微分方程的数值解法是十分必要的。,用Matlab软件求常微分方程的数值解,t，x=solver（f,ts,x0,options）,1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.,2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组
21、.,注意:,例6-10 设有初值问题，试求其数值解，并与精确解相比较。(1)建立函数文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2)求解微分方程。t0=0;tf=10;y0=2;t,y=ode23(funt,t0,tf,y0);%求数值解y1=sqrt(t+1)+1;%求精确解tyy1,例6-11 求解著名的Van der Pol方程,并绘制系统时间响应曲线与系统相平面曲线。,首先将高阶微分方程必须等价地变换成一阶微分方程组.选择状态变量：,得方程组：,建立函数文件vdpol.m如下,function xdot=vdpol(t,x)xd
22、ot(1)=-(x(2)2-1)*x(1)-x(2);xdot(2)=x(1);xdot=xdot;,求解微分方程，命令如下：t0=0;tf=20;x0=0;0.25;t,x=ode45(vdpol,t0,tf,x0);t,x,结果为：ans=0 0 0.2500 0.0002-0.0001 0.2500 0.0004-0.0001 0.2500 0.0006-0.0002 0.2500.19.9121-0.7520 1.5581 20.0000-0.7961 1.4900,第一列为t的采样点，第2列与第3列分别为x与x对应于t点的值,做系统时间响应曲线与系统相平面曲线。subplot(1,2
23、,1);plot(t,x);title(系统时间响应曲线);xlabel(t);ylabel(x);subplot(1,2,2);plot(x(:,1),x(:,2);title(系统相平面曲线);xlabel(x(1);ylabel(x(2);,例6-12 有Lorenz模型的状态方程如下，试绘制系统时间响应曲线与系统相平面图。,其中,建立Lorenz模型的函数文件lorenz.m如下,function xdot=lorenz(t,x)xdot(1)=-8/3*x(1)+x(2)*x(3);xdot(2)=-10*x(2)+10*x(3);xdot(3)=-x(1)*x(2)+28*x(2)
24、-x(3);xdot=xdot;,求解微分方程，命令如下：x0=0,0,eps;t,x=ode23(lorenz,0,100,x0);t,x,结果为：0 0 0 0.0000 10.0000-0.0000 0.0000-0.0000 10.6250-0.0000-0.0000 0.0000 10.7813-0.0000 0.0000-0.0000.99.9877 29.4980 4.3925-0.9899 100.0000 28.4974 3.7674-1.0279,第一列为t的采样点，第2、3、4列分别为x1、x2与x3对应于t点的值,做系统时间响应曲线与系统相平面曲线。subplot(1,
25、2,1);plot(t,x);title(系统时间响应曲线);xlabel(t);ylabel(x);subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3);title(系统相平面曲线);xlabel(x(1);ylabel(x(2);zlabel(x(3);,引例,标准形式：,求解无约束最优化问题的基本思想,求解的基本思想(以二元函数为例),5,3,1,连续可微,多局部极小,唯一极小(全局极小),2.优化函数的输入变量,使用优化函数或优化工具箱中其它优化函数时,输入变量见下表:,3.优化函数的输出变量下表:,6.4 最优化问题求解 6.4.1 无约束最优化问题求解
26、MATLAB提供了3个求最小值的函数，它们的调用格式为：（1）x,fval=fminbnd(fname,x1,x2,options)：求一元函数在(xl，x2)区间中的极小值点x和最小值fval。（2）x,fval=fminsearch(fname,x0,options)：基于单纯形算法求多元函数的极小值点x和最小值fval。（3）x,fval=fminunc(fname,x0,options)：基于拟牛顿法求多元函数的极小值点x和最小值fval。,控制参数options的设置,(3)MaxIter:允许进行迭代的最大次数,取值为正整数.,Options中常用的几个参数的名称、含义、取值如下:
27、,(1)Display:显示水平.取值为off时,不显示输出;取值为iter时,显示每次迭代的信息;取值为final时,显示最终结果.默认值为final.,(2)MaxFunEvals:允许进行函数评价的最大次数,取值为正整数.,例：opts=optimset(Display,iter,TolFun,1e-8)该语句创建一个称为opts的优化选项结构,其中显示参数设为iter,TolFun参数设为1e-8.,控制参数options可以通过函数optimset创建或修改。命令的格式如下：,(1)options=optimset(optimfun)创建一个含有所有参数名,并与优化函数optimfu
28、n相关的默认值的选项结构options.,（2）options=optimset(param1,value1,param2,value2,.)创建一个名称为options的优化选项参数,其中指定的参数具有指定值,所有未指定的参数取默认值.,(3)options=optimset(oldops,param1,value1,param2,value2,.)创建名称为oldops的参数的拷贝,用指定的参数值修改oldops中相应的参数.,返回,用Matlab解无约束优化问题,其中（3）、（4）、（5）的等式右边可选用（1）或（2）的等式右边。函数fminbnd的算法基于黄金分割法和二次插值法，它要求
29、目标函数必须是连续函数，并可能只给出局部最优解。,常用格式如下：（1）x=fminbnd(fun,x1,x2)（2）x=fminbnd(fun,x1,x2，options)（3）x，fval=fminbnd（.）（4）x，fval，exitflag=fminbnd（.）（5）x，fval，exitflag，output=fminbnd（.）,例6-13 求f(x)=x3-2x-5在0,5内的最小值点。(1)建立函数文件mymin.m。function fx=mymin(x)fx=x.3-2*x-5;(2)调用fmin函数求最小值点。x=fminbnd(mymin,0,5)x=0.8165,例6
30、-13 求f(x)=x3-2x-5在0,5内的最小值点和最大值点(1)建立M文件如下mymin.m。f=x3-2*x-5;%注意函数要命名为字符串fplot(f,0,5);%在某个区间作一元函数图语句xmin,ymin=fminbnd(f,0,5)f1=-x3+2*x+5;xmax,ymax=fminbnd(f1,0,5)xmin=0.8165 ymin=-6.0887xmax=5.0000ymax=109.9967%MATLAB实际算出的是-109.9967,命令格式为:（1）x=fminunc（fun,X0）；或x=fminsearch（fun,X0）（2）x=fminunc（fun,X0
31、，options）；或x=fminsearch（fun,X0，options）（3）x，fval=fminunc（.）；或x，fval=fminsearch（.）（4）x，fval，exitflag=fminunc（.）；或x，fval，exitflag=fminsearch（5）x，fval，exitflag，output=fminunc（.）；或x，fval，exitflag，output=fminsearch（.）,2、多元函数无约束优化问题,标准型为：min F(X),3 fminunc为中型优化算法的步长一维搜索提供了两种算法，由options中参数LineSearchType控制：
32、LineSearchType=quadcubic(缺省值)，混合的二次和三 次多项式插值；LineSearchType=cubicpoly，三次多项式插,使用fminunc和 fminsearch可能会得到局部最优解.,说明:,fminsearch是用单纯形法寻优.fminunc的算法见以下几点说明：,1 fminunc为无约束优化提供了大型优化和中型优化算法。由options中的参数LargeScale控制：LargeScale=on(默认值),使用大型算法LargeScale=off(默认值),使用中型算法,2 fminunc为中型优化算法的搜索方向提供了4种算法，由 options中的参
33、数HessUpdate控制：HessUpdate=bfgs（默认值），拟牛顿法的BFGS公式；HessUpdate=dfp，拟牛顿法的DFP公式；HessUpdate=steepdesc，最速下降法,例3 min f(x)=(4x12+2x22+4x1x2+2x2+1)*exp(x1),1、编写M-文件 fun1.m:function f=fun1(x)f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);2、输入M文件wliti3.m如下:x0=-1,1;x=fminunc(fun1,x0);y=fun1(x),3、运行结果:x=0.5000-1
34、.0000 y=1.3029e-10,1.为获得直观认识，先画出Rosenbrock 函数的三维图形,输入以下命令：x,y=meshgrid(-2:0.1:2,-1:0.1:3);z=100*(y-x.2).2+(1-x).2;mesh(x,y,z),2.画出Rosenbrock 函数的等高线图,输入命令：contour(x,y,z,20)hold on plot(-1.2,2,o);text(-1.2,2,start point)plot(1,1,o)text(1,1,solution),Rosenbrock 函数 f（x1，x2）=100（x2-x122+(1-x1)2 的最优解（极小）为
35、x*=（1，1），极小值为f*=0.试用 不同算法（搜索方向和步长搜索）求数值最优解.初值选为x0=（-1.2,2）.,3.用fminsearch函数求解,输入命令:f=100*(x(2)-x(1)2)2+(1-x(1)2;x,fval,exitflag,output=fminsearch(f,-1.2 2),运行结果:x=1.0000 1.0000fval=1.9151e-010exitflag=1output=iterations:108 funcCount:202 algorithm:Nelder-Mead simplex direct search,4.用fminunc 函数,(1)建
36、立M-文件fun2.m function f=fun2(x)f=100*(x(2)-x(1)2)2+(1-x(1)2,(2)主程序wliti44.m oldoptions=optimset(fminunc)options=optimset(oldoptions,LargeScale,off)options11=optimset(options,HessUpdate,dfp)x11,fval11,exitflag11,output11=fminunc(fun2,-1.2 2,options11)pause options12=optimset(options,HessUpdate,dfp,Lin
37、eSearchType,cubicpoly)x12,fval12,exitflag12,output12=fminunc(fun2,-1.2 2,options12)pause,(2)主程序wliti44.m,options21=optimset(options,HessUpdate,bfgs)x21,fval21,exitflag21,output21=fminunc(fun2,-1.2 2,options21)pause options22=optimset(options,HessUpdate,bfgs,LineSearchType,cubicpoly)x22,fval22,exitfl
38、ag22,output22=fminunc(fun2,-1.2 2,options22)pause options31=optimset(options,HessUpdate,steepdesc)x31,fval31,exitflag31,output31=fminunc(fun2,-1.2 2,options31)pause options32=optimset(options,HessUpdate,steepdesc,MaxIter,8000,MaxFunEvals,8000)x32,fval32,exitflag32,output32=fminunc(fun2,-1.2 2,option
39、s32)pause options33=optimset(options,HessUpdate,steepdesc,MaxIter,9000,MaxFunEvals,9000)x33,fval33,exitflag33,output33=fminunc(fun2,-1.2 2,options33),Rosenbrock函数不同算法的计算结果,可以看出，最速下降法的结果最差.因为最速下降法特别不适合于从一狭长通道到达最优解的情况.,6.4.2 有约束最优化问题求解MATLAB最优化工具箱提供了一个fmincon函数，专门用于求解各种约束下的最优化问题。该函数的调用格式为：x,fval=fminc
40、on(fname,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,options)其中，x、fval、fname、x0和options的含义与求最小值函数相同。其余参数为约束条件，参数NonF为非线性约束函数的M文件名。如果某个约束不存在，则用空矩阵来表示。,用MATLAB软件求解,其输入格式如下:1.x=quadprog(H,C,A,b);2.x=quadprog(H,C,A,b,Aeq,beq);3.x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB);4.x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,X0);5.x=quadprog(
41、H,C,A,b,Aeq,beq,VLB,VUB,X0,options);6.x,fval=quaprog(.);7.x,fval,exitflag=quaprog(.);8.x,fval,exitflag,output=quaprog(.);,1、二次规划,1.首先建立M文件fun.m,定义目标函数F（X）:function f=fun(X);f=F(X);,2、一般非线性规划,其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同.用Matlab求解上述问题，基本步骤分三步：,3.建立主程序.非线性规划求解的函数是fmincon,命令的
42、基本格式如下：(1)x=fmincon(fun,X0,A,b)(2)x=fmincon(fun,X0,A,b,Aeq,beq)(3)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB)(4)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon)(5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options)(6)x,fval=fmincon(.)(7)x,fval,exitflag=fmincon(.)(8)x,fval,exitflag,output=fmincon(.),输出极值点,M
43、文件,迭代的初值,参数说明,变量上下限,注意：1 fmincon函数提供了大型优化算法和中型优化算法。默认时，若在fun函数中提供了梯度（options参数的GradObj设置为on），并且只有上下界存在或只有等式约束，fmincon函数将选择大型算法。当既有等式约束又有梯度约束时，使用中型算法。2 fmincon函数的中型算法使用的是序列二次规划法。在每一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日Hessian矩阵。3 fmincon函数可能会给出局部最优解，这与初值X0的选取有关。,1、写成标准形式：s.t.,2x1+3x2 6 s.t x1+4x2 5 x1,x2 0,例2,
44、2、先建立M-文件 fun3.m:function f=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)2,3、再建立主程序youh2.m：x0=1;1;A=2 3;1 4;b=6;5;Aeq=;beq=;VLB=0;0;VUB=;x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB),4、运算结果为：x=0.7647 1.0588 fval=-2.0294,1先建立M文件 fun4.m,定义目标函数:function f=fun4(x);f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2
45、*x(2)+1);,x1+x2=0 s.t.1.5+x1x2-x1-x2 0-x1x2 10 0,例3,2再建立M文件mycon.m定义非线性约束：function g,ceq=mycon(x)g=x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;,3主程序youh3.m为:x0=-1;1;A=;b=;Aeq=1 1;beq=0;vlb=;vub=;x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon),3.运算结果为：x=-1.2250 1.2250 fval=1.8951,例4,1先建立M-文件fun
46、.m定义目标函数:function f=fun(x);f=-2*x(1)-x(2);,2再建立M文件mycon2.m定义非线性约束：function g,ceq=mycon2(x)g=x(1)2+x(2)2-25;x(1)2-x(2)2-7;,3.主程序fxx.m为:x0=3;2.5;VLB=0 0;VUB=5 10;x,fval,exitflag,output=fmincon(fun,x0,VLB,VUB,mycon2),4.运算结果为:x=4.0000 3.0000fval=-11.0000exitflag=1output=iterations:4 funcCount:17 stepsiz
47、e:1 algorithm:1x44 char firstorderopt:cgiterations:,6.4.3 线性规划问题求解在MATLAB中求解线性规划问题使用函数linprog，其调用格式为x,fval=linprog(f,A,b,Aeq,beq,lbnd,ubnd)其中，x是最优解，fval是目标函数的最优值。函数中的各项参数是线性规划问题标准形式中的对应项，x、b、beq、lb、ub是向量，A、Aeq为矩阵，f为目标函数系数向量。,用MATLAB优化工具箱解线性规划,命令：x=linprog（c，A，b）,2、模型：min z=cX,命令：x=linprog（c，A，b，Aeq,
48、beq）,注意：若没有不等式：存在，则令A=，b=.,命令：1 x=linprog（c，A，b，Aeq,beq,VLB，VUB）2 x=linprog（c，A，b，Aeq,beq,VLB，VUB,X0）,注意：1 若没有等式约束:,则令Aeq=,beq=.2其中X0表示初始点,4、命令：x,fval=linprog()返回最优解及处的目标函数值fval.,解 编写M文件xxgh1.m如下：c=-0.4-0.28-0.32-0.72-0.64-0.6;A=0.01 0.01 0.01 0.03 0.03 0.03;0.02 0 0 0.05 0 0;0 0.02 0 0 0.05 0;0 0 0.03 0 0 0.08;b=850;700;100;900;Aeq=;beq=;vlb=0;0;0;0;0;0;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),解:编写M文件xxgh2.m如下：c=6 3 4;A=0 1 0;b=50;Aeq=1 1 1;beq=120;vlb=30,0,20;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),
