第九章最优化方法

1、第九章 最优化方法 本章主要介绍线性规划、0-1规划、非线性规划等问题的MATLAB求解。 9.1 线性规划（Linear Programming，简写为LP）问题线性规划问题就是求多变量线性函数在线性约束条件下的最优值。满足约束条件的解称为可行解，所有可行解构成的集合称为可行域，满足目标式的可行解称为最优解。MATLAB解决的线性规划问题的标准形式为：min 其中为列向量，为矩阵。其它形式的线性规划问题都可经过适当变换化为此标准形式。在MATLAB中求解线性规划问题函数为linprog，其使用格式为：x, fval, exitflag, output, lambda = linprog(f,

2、 A, b, Aeq, beq, lb, ub)输入部分：其中各符号对应线性规划问题标准形式中的向量和矩阵，如果约束条件中有缺少，则其相应位置用空矩阵代替。输出部分：其中x为最优解，用列向量表示；fval为最优值；exitflag为退出标志，若exitflag=1表示函数有最优解，若exitflag=0表示超过设定的迭代最大次数，若exitflag=-2，表示约束区域不可行，若exitflag=-3，表示问题无解，若exitflag=-4，表示执行迭代算法时遇到NaN，若exitflag=-5，表示原问题和对偶问题均不可行，若exitflag=-7，表示搜索方向太小，不能继续前进；output

3、表明算法和迭代情况；lambda表示存储情况。例1 用linprog函数求下面的线性规划问题输入如下：>> f = -5, -4,-6;A = 1 -1 1; 3 2 4; 3 2 0;b = 20; 42; 30;lb = zeros(3,1);x,fval,exitflag,output,lambda = linprog(f,A,b,lb)注意：由于该问题没有等式约束，所以输入格式中相应的位置用代替，变量没有上限约束，所以ub也用代替，但由于其在最后，可以不写。输出结果如下：Optimization terminated.x = % 最优解 0.0000 15.0000 3.0

4、000fval = %最优值 -78.0000exitflag = %函数收敛于解 1output = iterations: 6 algorithm: 'large-scale: interior point' cgiterations: 0 message: 'Optimization terminated.'lambda = ineqlin: 3x1 double eqlin: 0x1 double upper: 3x1 double lower: 3x1 double例2 一家家具公司生产桌子和椅子，用于生产的全部劳动力共计450个工时，原料是400个单

5、位的木材。每张桌子使用15个工时的劳动力，20个单位的木材，售价为80元。每张椅子使用10个工时，用料5个单位，售价45元。问为达到最大效益，应如何安排生产？解 设生产桌子个，椅子个，建立如下模型：输入如下：>> f = -80,-45;A = 20, 5; 15, 10;b = 400, 450;lb = zeros(2,1);x, fval, exitflag = linprog(f,A,b,lb)结果如下：Optimization terminated.x = 14.0000 24.0000fval = -2.2000e+003exitflag = 1注意：由于linprog

6、是求目标函数的最小值，如求目标函数的最大值，可先求出的最小值fval，则-fval就是的最大值。本例只输出最优解、最优值和退出标志，可见生产14个桌子，24个椅子，可获得最大利润2200元。9.2 0-1规划0-1规划是一种特殊形式的整数规划，它的决策变量仅取值0或1.一般用0表示放弃，1表示选取，故0-1规划可以用来处理选址问题、指派问题、装箱问题、项目评价、资金分配、生产计划安排等问题。在MATALB中求解0-1规划问题函数为bintprog，其针对下述0-1规划：min 其中为列向量，为矩阵。使用格式为：x, fval, exitflag, output = bintprog(f, A,

7、 b, Aeq, beq)输入部分：其中各符号对应0-1规划问题标准形式中的向量和矩阵，如果约束条件中有缺少，则其相应位置用空矩阵代替。输出部分：其中x为最优解，用列向量表示；fval为最优值；exitflag为退出标志，若exitflag=1表示函数有最优解，若exitflag=-1表示超过设定的迭代最大次数，若exitflag=-2，表示问题不可行，若exitflag=-4，表示搜索节点数超过了设定的搜索节点最大个数，若exitflag=-5，表示搜索时间超过了设定的指令运行的最大秒数，若exitflag=-6，表示LP求解器在某节点处求解LP松弛问题时的迭代次数超过了设定的迭代次数；ou

8、tput包含使用算法、迭代次数、搜索过的节点数、算法执行时间、算法终止原因。例3 求解下述0-1规划问题。利用bintprog函数求解如下：>> f=-1;2;2;-6;-4;A=3,2,-1,4,2; 2,4, -2,-1,-2;b=5;5 ;x,fval,exit,output=bintprog(f,A,b)输出结果：Optimization terminated.x = 1 1 1 0 0fval = -5exit = 1output = iterations: 4 nodes: 1 time: 0.0781 algorithm: 'LP-based branch-a

9、nd-bound' branchStrategy: 'maximum integer infeasibility' nodeSrchStrategy: 'best node search' message: 'Optimization terminated.这表明，当时，目标函数取最大值5。例4 资金分配问题，某企业在今后3年内有5个工程考虑开工，每项工程的期望收入和年度费用如表所示。假定每一项已经批准的工程要在整3年内完成。企业应如何选择工程，使企业总收入最大？工程费用（千元）收入（千元）第一年第二年第三年1518202471040339220

10、4741155861030最大可用资金数252525解 设决策变量为：其中，表示放弃第项工程，表示选择第项工程。该问题的数学模型为：利用bintprog函数求解如下：f=-20;40;20;15;30;A=5,4,3,7,8;1,7,9,4,6;8,10,2,1,10;b=25;25;25;x,fv,exit=bintprog(f,A,b)输出结果如下：Optimization terminated.x = 1 1 1 1 0fv = -95exit = 1上述结果表明，该企业选择第1，第2，第3，第4项工程，可以获得最大利润95千元。指派问题：设有项工作分配给个人去做，每人做一项，由于个人的

11、工作效率不同，因而完成同一工作所需时间也不同，设人员完成工作所需时间为（称为效率矩阵），问如何分配工作，使得完成所有工作所用的总时间最少？这类问题称为指派问题（Assignment Problem），也称最优匹配问题，它是一类典型的0-1规划问题。例5 有4个工人被分派做4项工作，每项工作只能一人做，每人只能做一项工作。现设每个工人做每项工作的耗时如表所示，求总耗时最少的分配方案。表1时间表 时间单位：小时 工作耗时工人1234115182124219232218326171619419212317解 设决策变量为，只取0或1。表示工人完成工作；表示工人不做工作。于是，上述问题的数学模型如下：

12、下面给出Matlab计算程序。>> f=15;18;21;24;19;23;22;18;26;17;16;19;19;21;23;17;o4=ones(1,4);z4=zeros(1,4);z8=zeros(1,8);z12=zeros(1,12);Aeq1=o4,z12;z4,o4,z8;z8,o4,z4;z12,o4;e4=eye(4);Aeq2=e4,e4,e4,e4;Aeq = Aeq1;Aeq2;beq=ones(8,1);x,fv,exit=bintprog(f,Aeq,beq);x=transpose(reshape(x,4,4)计算接可以得出：fv = 70exit

13、 = 1指派方阵x = 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0它表示的最优指派方案为：工人1做工作1；工人2做工作4；工人3做工作3；工人4做工作2。该指派方案耗时最少，为70小时。9.3 非线性规划9.3.1 有约束的一元函数的极（最）小值在MATLAB中使用fminbnd函数求一元函数在区间上的极（最）小值。其调用格式有：x,fval,exitflag,output = fminbnd(fun,x1,x2,options)其中输入参数fun为目标函数的表达式字符串或MATLAB定义的M函数，x1和x2为自变量x的取值范围，options为指定优化参数选项；输出参数x

14、为函数fun在区间上取最小值时的x值，fval为目标函数的最小值，exitflag为终止迭代的条件，output为优化信息，参见optimset函数。若参数exitflag>0，表示函数收敛于x，若exitflag=0，表示超过函数估计值或迭代的最大数字，exitflag<0表示函数不收敛于x；若参数output=iterations表示迭代次数，output=funccount表示函数赋值次数，output=algorithm表示所使用的算法。例6 计算函数在区间(0,1)内的最小值。在命令窗口输入：>> x,fval,exitflag,output=fminbnd(

15、'(x3+cos(x)+x*log(x)/exp(x)',0,1)输出结果为：x = 0.5223fval = 0.3974exitflag = 1output = iterations: 8 funcCount: 9 algorithm: 'golden section search, parabolic interpolation' message: 1x112 char表明该函数在x =0.5223时取最小值0.3974。ezplot('(x3+cos(x)+x*log(x)/exp(x)',0,1)为简化输出结果，可输入：>>

16、 x,fval=fminbnd('(x3+cos(x)+x*log(x)/exp(x)',0,1)9.3.2 无约束多元函数极（最）小值在MATLAB中使用fminsearch函数和fminunc函数求无约束多元函数的极（最）小值。（1）fminsearch函数该函数的调用格式为：x,fval,exitflag,output= fminsearch(fun,x0,options)其中输入参数fun为目标函数的表达式字符串或MATLAB自定义的M函数，x0为初始点，options为指定优化参数选项；输出参数x为多元函数的最小值点，fval为最小值，exitflag和output与

17、单变量情形一致。例7 求的最小值点。在命令窗口输入：>> x,fval=fminsearch('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2', 0,0)输出结果为：x = 1.0016 0.8335fval = -3.3241表明当x1=1.0016，x2=0.8335时函数的最小值为-3.3241。（2）fminunc函数该函数的调用格式为：x,fval,exitflag,output,grad,hessian = fminunc (fun,x0,options)其中输入参数fun为目标函数的表达式字符串或MATLAB自定义的M

18、函数，x0为初始点，options为指定优化参数选项；输出参数x为多元函数的最小值点，fval为最小值，exitflag和output与单变量情形一致，grad为函数在解x处的梯度值，hessian为目标函数在解x处的海赛（Hessian）值。fminsearch利用了单纯形法的原理，fminunc利用了拟牛顿法的原理。当函数的阶数大于2时，使用fminunc比fminsearch更有效，但当所选函数高度不连续时，使用fminsearch效果较好，这两个函数都容易陷入局部优化，并且结果的正确与否还要取决于初值点x0的选取。利用fiminunc求例7的最小值。在名利窗口输入：>>x,

19、fval=fminunc('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2',0.5,0.5)输出结果为：x = 1.0016 0.8335fval = -3.3241此结果和fminsearch函数的求解结果一样。9.3.3有约束的多元函数最小值非线性有约束的多元函数的标准形式为：其中：x,b,beq,lb,ub是向量，A,Aeq为矩阵，C(x)、Ceq(x)是返回向量的函数，f(x)为目标函数，f(x)、C(x)、Ceq(x)可以是非线性函数。在MATLAB用fmincon函数，求有约束的多元函数的最小值，其调用格式有：x = fmincon

20、(fun,x0,A,b)x = fmincon(fun,x0,A,b,Aeq,beq)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)x,fval = fmincon()x,fval,exitflag = fmincon()x,fval,exitflag,output = fmincon()x,fval,exitflag,output,lambda = fmincon()x,

21、fval,exitflag,output,lambda,grad = fmincon()x,fval,exitflag,output,lambda,grad,hessian = fmincon()其中fun为目标函数，它可用前面的方法定义；x0为初始值；A、b满足线性不等式约束，若没有不等式约束，则取A=，b=；Aeq、beq满足等式约束，若没有，则取Aeq=，beq=；lb、ub满足，若没有界，可设lb=，ub=；nonlcon的作用是通过接受的向量x来计算非线性不等约束和等式约束分别在x处的估计C和Ceq，通过M文件建立，如：>>x = fmincon(myfun,x0,A,b

22、,Aeq,beq,lb,ub,mycon)，先建立非线性约束函数，并保存为mycon.m：function C,Ceq = mycon(x)C = % 计算x处的非线性不等约束的函数值。Ceq = % 计算x处的非线性等式约束的函数值。输出参数x为多元函数的最小值点，fval为最小值，exitflag和output与单变量情形一致，lambda是Lagrange乘子，它体现哪一个约束有效，grad为函数在解x处的梯度值，hessian为目标函数在解x处的海赛（Hessian）值。例8 求下面问题在初始点x=(10, 10, 10)处的最优解。min sub.to 解：约束条件的标准形式为sub

23、.to 在命令窗口输入：>> fun= '-x(1)*x(2)*x(3)'x0=10,10,10;A=-1 -2 -2;1 2 2;b=0;72;x,fval=fmincon(fun,x0,A,b)输出结果为：Optimization terminated: magnitude of directional derivative in search direction less than 2*options.TolFun and maximum constraint violation is less than options.TolCon.Active inequ

24、alities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 2 x = 24.0000 12.0000 12.0000fval = -3.4560e+003 例9 求表面积为36而体积最大的长方体的体积。解 设长方体的三个棱长分别为x1,x2,x3，则体积v=x1*x2*x3，且满足约束条件2(x1*x2+x2*x3+x3*x1)=36。首先建立目标函数的M文件myopt1.m，内容如下：function f = myopt1(x)f =-x(1)*x(2)*x(3); % 负号表示最小值，即x(1)*

25、x(2)*x(3)的最大值再建立非线性约束函数的M文件mycon.m，内容如下：function C,Ceq = mycon(x)C = ;Ceq = x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-18;由于没有非线性不等式，故C为。在命令窗口输入：>> x,fval = fmincon('myopt1',1 1 1,0 0 0,'mycon')输出结果为：x = 2.4495 2.4495 2.4495fval = -14.6969结果表明当x1=x2=x3=2.4495时，体积最大值为14.6969。该问题的精确解为x1=x2=x3

26、=时最大体积为。9.4 二次规划二次规划问题（quadratic programming）的标准形式为：sub.to 其中，H、A、Aeq为矩阵，f、b、beq、lb、ub、x为向量。其它形式的二次规划问题都可转化为标准形式。MATLAB中使用quadprog函数求解二次规划，其调用格式有：x = quadprog(H,f,A,b) x = quadprog(H,f,A,b,Aeq,beq) x = quadprog(H,f,A,b,Aeq,beq,lb,ub) x = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)x = quadprog(H,f,A,b,Aeq,beq

27、,lb,ub,x0,options)x,fval = quadprog() x,fval,exitflag = quadprog() x,fval,exitflag,output = quadprog()x,fval,exitflag,output,lambda = quadprog()例10 求解下面二次规划问题解：则，在命令窗口中输入：>> H = 1 -1; -1 2 ;f = -2; -6;A = 1 1; -1 2; 2 1;b = 2; 2; 3;lb = zeros(2,1);x,fval,exitflag,output,lambda = quadprog(H,f,A

28、,b, , ,lb)输出结果为：x = %最优解 0.6667 1.3333fval = %最优值 -8.2222exitflag = %收敛 1output = iterations: 3 algorithm: 'medium-scale: active-set' firstorderopt: cgiterations: lambda = lower: 2x1 double upper: 2x1 double eqlin: 0x1 double ineqlin: 3x1 double>> lambda.ineqlinans = 3.1111 0.4444 0>

29、;> lambda.lowerans = 0 0此结果说明第1、2个约束条件有效，其余无效。9.5 多目标规划多目标规划是指在一组约束下，对多个不同目标函数进行优化。它的一般形式为其中：x、b、beq、lb、ub是向量；A、Aeq为矩阵；C(x)、Ceq(x)和F(x)是返回向量的函数；F(x)、C(x)、Ceq(x)可以是非线性函数；weight为权值系数向量，用于控制对应的目标函数与用户定义的目标函数值的接近程度；goal为用户设计的与目标函数相应的目标函数值向量；为一个松弛因子标量；F(x)为多目标规划中的目标函数向量。MATLAB中求解多目标规划的函数为fgoalattain，其

30、调用格式有：x = fgoalattain(fun,x0,goal,weight)x = fgoalattain(fun,x0,goal,weight,A,b)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options)x,fval = fgoalattain()x,fval,attainfactor = fgoalattain()x,fval,attainfactor,exitflag = fgoalattain()x,fval,attainfactor,exitflag,output = fgoalattain()x,fval,attainfactor,exitflag,output,lambda = fgoalattain()输入参数x0为初始解向量；fun为多目标函数的文件名字符