matlab处理优化问题

176
第5章优化问题
5.1线性规划问题
线性规划问题是目标函数和约束条件均为线性函数的问题，MATLAB6.0解决的线性规划问题的标准形式为：min n
R x x f ∈'sub.to ：b x A ≤⋅beq x Aeq =⋅ub
x lb ≤≤其中f 、x 、b 、beq 、lb 、ub 为向量，A 、Aeq 为矩阵。

其它形式的线性规划问题都可经过适当变换化为此标准形式。

在MATLAB6.0版中，线性规划问题（Linear Programming ）已用函数linprog 取代了MATLAB5.x 版中的lp 函数。

当然，由于版本的向下兼容性，一般说来，低版本中的函数在6.0版中仍可使用。

函数linprog
格式x =linprog(f,A,b)%求min f '*x sub.to b x A ≤⋅线性规划的最优解。

x =linprog(f,A,b,Aeq,beq)%等式约束beq x Aeq =⋅，若没有不等式约
束b x A ≤⋅，则A=[]，b=[]。

x =linprog(f,A,b,Aeq,beq,lb,ub)%指定x 的范围ub x lb ≤≤，若没有等
式约束beq x Aeq =⋅，则Aeq=[]，beq=[]
x =linprog(f,A,b,Aeq,beq,lb,ub,x0)%设置初值x0x =linprog(f,A,b,Aeq,beq,lb,ub,x0,options)%options 为指定的优化
参数
[x,fval]=linprog(…)%返回目标函数最优值，即fval=f '*x 。

[x,lambda,exitflag]=linprog(…)%lambda 为解x 的Lagrange 乘子。

[x,lambda,fval,exitflag]=linprog(…)%exitflag 为终止迭代的错误
条件。

[x,fval,lambda,exitflag,output]=linprog(…)%output 为关于优化的
一些信息
说明若exitflag>0表示函数收敛于解x ，exitflag=0表示超过函数估值或迭代的最大数字，exitflag<0表示函数不收敛于解x ；若lambda=lower 表示下界lb ，lambda=upper 表示上界ub ，lambda=ineqlin 表示不等式约束，lambda=eqlin 表示
177
等式约束，lambda 中的非0元素表示对应的约束是有效约束；output=iterations 表示迭代次数，output=algorithm 表示使用的运算规则，output=cgiterations 表示PCG 迭代次数。

例5-1求下面的优化问题
min 321x 6x 4x 5---sub.to 20
x x x 321≤+-42
x 4x 2x 3321≤++30
x 2x 321≤+3
21x 0,x 0,x 0≤≤≤解：>>f =
[-5;-4;-6];>>A =[1-11;324;3
20];
>>b =[20;42;30];
>>lb =zeros(3,1);>>[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb)
结果为：x =
%最优解
0.000015.00003.0000fval =%
最优值
-78.0000exitflag =%
收敛
1output =
iterations:6%迭代次数
cgiterations:0algorithm:'lipsol'%所使用规则
lambda =
ineqlin:[3x1double]
eqlin:[0x1double]
upper:[3x1double]lower:[3x1double]
>>lambda.ineqlin ans =0.00001.5000
0.5000
>>lambda.lower ans =1.00000.00000.0000
表明：不等约束条件2和3以及第1个下界是有效的
178
5.2foptions 函数
对于优化控制，MATLAB 提供了18个参数，这些参数的具体意义为：
options(1)-参数显示控制（默认值为0）。

等于1时显示一些结果。

options(2)-优化点x 的精度控制(默认值为1e-4)。

options(3)-优化函数F 的精度控制(默认值为1e-4)。

options(4)-违反约束的结束标准(默认值为1e-6)。

options(5)-算法选择，不常用。

options(6)-优化程序方法选择，为0则为BFCG 算法，为1则采用DFP
算法。

options(7)-线性插值算法选择，为0则为混合插值算法，为1则采用立
方插算法。

options(8)-函数值显示(目标—达到问题中的Lambda )options(9)-若需要检测用户提供的梯度，则设为1。

options(10)-函数和约束估值的数目。

options(11)-函数梯度估值的个数。

options(12)-约束估值的数目。

options(13)-等约束条件的个数。

options(14)-函数估值的最大次数（默认值是100×变量个数）options(15)-用于目标—达到问题中的特殊目标。

options(16)-优化过程中变量的最小有限差分梯度值。

options(17)-优化过程中变量的最大有限差分梯度值。

options(18)-步长设置(默认为1或更小)。

Foptions 已经被optimset 和optimget 代替，详情请查函数optimset 和optimget 。

5.3
非线性规划问题5.3.1有约束的一元函数的最小值
单变量函数求最小值的标准形式为)x (f min x sub.to 2
1x x x <<在MATLAB5.x 中使用fmin 函数求其最小值。

函数fminbnd
格式x =fminbnd(fun,x1,x2)%返回自变量x 在区间21x x x <<上函数fun
取最小值时x 值，fun 为目标函数的表达
式字符串或MATLAB 自定义函数的函数
柄。

179
x =fminbnd(fun,x1,x2,options)%options 为指定优化参数选项[x,fval]=fminbnd(…)%fval 为目标函数的最小值[x,fval,exitflag]=fminbnd(…)%xitflag 为终止迭代的条件[x,fval,exitflag,output]=fminbnd(…)%output 为优化信息
说明若参数exitflag>0，表示函数收敛于x ，若exitflag=0，表示超过函数估计值或迭代的最大数字，exitflag<0表示函数不收敛于x ；若参数output=iterations 表示迭代次数，output=funccount 表示函数赋值次数，output=algorithm 表示所使用的算法。

例5-2计算下面函数在区间(0，1)内的最小值。

x
3e x log x x cos x )x (f ++=解：>>[x,fval,exitflag,output]=fminbnd('(x^3+cos(x)+x*log(x))/exp(x)',0,1)x =0.5223fval =0.3974exitflag =1output =iterations:9funcCount:9algorithm:'golden section search,parabolic interpolation'
例5-3在[0，5]上求下面函数的最小值1
)3x ()x (f 3--=解：先自定义函数：在MATLAB 编辑器中建立M 文件为：function f =myfun(x)f =(x-3).^2-1;
保存为myfun.m ，然后在命令窗口键入命令：>>x=fminbnd(@myfun,0,5)
则结果显示为：x =3
5.3.2无约束多元函数最小值
多元函数最小值的标准形式为)x (f min x
其中：x 为向量，如]
x ,,x ,x [x n 21 =在MATLAB5.x 中使用fmins 求其最小值。

命令利用函数fminsearch 求无约束多元函数最小值
函数fminsearch
格式x =fminsearch(fun,x0)%x0为初始点，fun 为目标函数的表达式字符
串或MATLAB 自定义函数的函数柄。

x =fminsearch(fun,x0,options)%options 查optimset
180
[x,fval]=fminsearch(…)%最优点的函数值[x,fval,exitflag]=fminsearch(…)%exitflag 与单变量情形一致[x,fval,exitflag,output]=fminsearch(…)%output 与单变量情形一致
注意：fminsearch 采用了Nelder-Mead 型简单搜寻法。

例5-4求222132131x x x 10x x 4x 2y +-+=的最小值点
解：>>X=fminsearch('2*x(1)^3+4*x(1)*x(2)^3-10*x(1)*x(2)+x(2)^2',[0,0])结果为X =1.00160.8335
或在MATLAB 编辑器中建立函数文件function f=myfun(x)f=2*x(1)^3+4*x(1)*x(2)^3-10*x(1)*x(2)+x(2)^2;
保存为myfun.m ，在命令窗口键入>>X=fminsearch ('myfun',[0,0])或>>X=fminsearch(@myfun,[0,0])结果为：X =1.00160.8335
命令利用函数fminunc 求多变量无约束函数最小值
函数fminunc
格式x =fminunc(fun,x0)%返回给定初始点x0的最小函数值点
x =fminunc(fun,x0,options)%options 为指定优化参数[x,fval]=fminunc(…)%fval 最优点x 处的函数值[x,fval,exitflag]=fminunc(…)%exitflag 为终止迭代的条件，与上
同。

[x,fval,exitflag,output]=fminunc(…)%output 为输出优化信息[x,fval,exitflag,output,grad]=fminunc(…)%grad 为函数在解x 处的
梯度值
[x,fval,exitflag,output,grad,hessian]=fminunc(…)%目标函数在解x
处的海赛（Hessian ）
值
注意：当函数的阶数大于2时，使用fminunc 比fminsearch 更有效，但当所选函数高度不连续时，使用fminsearch 效果较好。

例5-5求222121x x x 2x 3)x (f ++=的最小值。

>>fun='3*x(1)^2+2*x(1)*x(2)+x(2)^2';>>x0=[11];>>[x,fval,exitflag,output,grad,hessian]=fminunc(fun,x0)
结果为：x =1.0e-008*-0.75910.2665fval =
181
1.3953e-016exitflag =1output =iterations:3funcCount:16stepsize:1.2353firstorderopt:1.6772e-007algorithm:'medium-scale:Quasi-Newton line search'grad =1.0e-006*-0.16770.0114hessian =6.0000
2.00002.0000 2.0000
或用下面方法：>>fun=inline('3*x(1)^2+2*x(1)*x(2)+x(2)^2')fun =Inline function:fun(x)=3*x(1)^2+2*x(1)*x(2)+x(2)^2>>x0=[11]；>>x=fminunc(fun,x0)x =1.0e-008*-0.75910.2665
5.3.3有约束的多元函数最小值
非线性有约束的多元函数的标准形式为：
)x (f min x
sub.to 0
)x (C ≤0
)x (Ceq =b
x A ≤⋅beq
x Aeq =⋅ub
x lb ≤≤其中：x 、b 、beq 、lb 、ub 是向量，A 、Aeq 为矩阵，C(x)、Ceq(x)是返回向量的函数，f(x)为目标函数，f(x)、C(x)、Ceq(x)可以是非线性函数。

在MATLAB5.x 中，它的求解由函数constr 实现。

函数fmincon
格式x =fmincon(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)
182
[x,fval]=fmincon(…)[x,fval,exitflag]=fmincon(…)[x,fval,exitflag,output]=fmincon(…)[x,fval,exitflag,output,lambda]=fmincon(…)[x,fval,exitflag,output,lambda,grad]=fmincon(…)[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(…)
参数说明：fun 为目标函数，它可用前面的方法定义；
x0为初始值；A 、b 满足线性不等式约束b x A ≤⋅，若没有不等式约束，则取A=[]，b=[]；
Aeq 、beq 满足等式约束beq x Aeq =⋅，若没有，则取Aeq=[]，
beq=[]；lb 、ub 满足ub x lb ≤≤，若没有界，可设lb=[]，ub=[]；
nonlcon 的作用是通过接受的向量x 来计算非线性不等约束0
)x (C ≤和等式约束0)x (Ceq =分别在x 处的估计C 和Ceq ，通过指
定函数柄来使用，如：>>x =
fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非
线性约束函数，并保存为mycon.m ：function [C,Ceq]=
mycon(x)
C =…
%计算x 处的非线性不等约束0)x (C ≤的函数值。

Ceq =…%计算x 处的非线性等式约束0)x (Ceq =的函数值。

lambda 是Lagrange 乘子，它体现哪一个约束有效。

output 输出优化信息；grad 表示目标函数在x 处的梯度；hessian 表示目标函数在x 处的Hessiab 值。

例5-6
求下面问题在初始点（0，1）处的最优解
min 2
1212221x 5x 2x x x x ---+sub.to 0
x )1x (221≥+--06x 3x 221≥+-解：约束条件的标准形式为
sub.to 0
x )1x (221≤--6
x 3x 221≤+-先在MATLAB 编辑器中建立非线性约束函数文件：function [c,ceq]=mycon (x)
183
c=(x(1)-1)^2-x(2);ceq=[];%无等式约束
然后，在命令窗口键入如下命令或建立M 文件：>>fun='x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2)';%目标函数>>x0=[01];>>A=[-23];%线性不等式约束>>b=6;>>Aeq=[];%无线性等式约束>>beq=[];>>lb=[];%x 没有下、上界>>ub=[];>>[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,@mycon)
则结果为x =34fval =-13exitflag =%解收敛1output =iterations:2funcCount:9stepsize:1algorithm:'medium-scale:SQP,Quasi-Newton,line-search'firstorderopt:[]cgiterations:[]lambda =lower:[2x1double]%x 下界有效情况，通过lambda.lower 可查看。

upper:[2x1double]%x 上界有效情况，为0表示约束无效。

eqlin:[0x1double]%线性等式约束有效情况，不为0表示约束有效。

eqnonlin:[0x1double]%非线性等式约束有效情况。

ineqlin:2.5081e-008%线性不等式约束有效情况。

ineqnonlin:6.1938e-008%非线性不等式约束有效情况。

grad =%目标函数在最小值点的梯度1.0e-006*-0.17760hessian =%目标函数在最小值点的Hessian 值1.0000-0.0000-0.0000 1.0000
例5-7求下面问题在初始点x=(10,10,10)处的最优解。

Min 3
21x x x )x (f -=Sub.to 72
x 2x 2x 0321≤++≤解：约束条件的标准形式为
sub.to
0x 2x 2x 321≤---72
x 2x 2x 321≤++>>fun='-x(1)*x(2)*x(3)';>>x0=[10,10,10];>>A=[-1-2-2;122];
184
>>b=[0;72];>>[x,fval]=fmincon(fun,x0,A,b)
结果为：x =24.000012.000012.0000fval =-3456
5.3.4二次规划问题
二次规划问题（quadratic programming ）的标准形式为：
x f x H x 2
1min '+'sub.to b x A ≤⋅beq
x Aeq =⋅b
u x b l ≤≤其中，H 、A 、Aeq 为矩阵，f 、b 、beq 、lb 、ub 、x 为向量
其它形式的二次规划问题都可转化为标准形式。

MATLAB5.x 版中的qp 函数已被6.0版中的函数quadprog 取代。

函数quadprog
格式x =quadprog(H,f,A,b)%其中H,f,A,b 为标准形中的参数，x 为目标
函数的最小值。

x =quadprog(H,f,A,b,Aeq,beq)%Aeq,beq 满足等约束条件
beq x Aeq =⋅。

x =quadprog(H,f,A,b,Aeq,beq,lb,ub)
%lb,ub 分别为解x 的下界与上
界。

x =quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)%x0为设置的初值x =quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)%options 为指定的
优化参数
[x,fval]=quadprog(…)%fval 为目标函数最优值[x,fval,exitflag]=quadprog(…)%exitflag 与线性规划中参数意义相
同
[x,fval,exitflag,output]=quadprog(…)
%output 与线性规划中参数意义相同
[x,fval,exitflag,output,lambda]=quadprog(…)
%lambda 与线性规划
中参数意义相同例5-8求解下面二次规划问题
21212221x 6x 2x x x x 21)x (f min ---+=sub.to
2
x x 21≤+
1852
x 2x 21≤+-3
x x 221≤+2
1x 0,x 0≤≤解：x
f x H x 21)x (f '+'=则⎥⎦⎤
⎢⎣⎡--=2111H ，⎥⎦⎤
⎢⎣⎡--=62f ，⎥⎦
⎤
⎢⎣⎡=21x x x 在MATLAB
中实现如下：>>H =[1-1;-12];
>>f =
[-2;-6];
>>A =[11;-12;21];
>>b =[2;2;3];>>lb =zeros(2,1);>>[x,fval,exitflag,output,lambda]=quadprog(H,f,A,b,[],[],lb)结果为：x =%最优解
0.66671.3333fval =%
最优值
-8.2222exitflag =%
收敛
1output =iterations:
3
algorithm:'medium-scale:active-set'
firstorderopt:[]cgiterations:[]
lambda =
lower:[2x1double]
upper:[2x1double]
eqlin:[0x1double]
ineqlin:[3x1double]>>lambda.ineqlin ans =3.11110.4444
>>lambda.lower ans =00
说明第1、2个约束条件有效，其余无效。

例5-9求二次规划的最优解
max f (x 1,x 2)=x 1x 2+3
sub.to x 1+x 2
-2=0解：化成标准形式：
3
x x )0,0(x x 0110)x x (213x x )x x (f min 2121212121-⎪⎭⎫
⎝⎛+⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛--=--=
186
sub.to x 1+x 2=2
在Matlab 中实现如下：>>H=[0,-1;-1,0];>>f=[0;0];>>Aeq=[11];>>b=2;>>[x,fval,exitflag,output,lambda]=quadprog(H,f,[],[],Aeq,b)
结果为：x = 1.00001.0000fval =-1.0000exitflag =1output =firstorderopt:0iterations:1cgiterations:1algorithm:[1x58char]lambda =eqlin:1.0000ineqlin:[]lower:[]upper:[]
5.4“半无限”有约束的多元函数最优解
“半无限”有约束多元函数最优解问题的标准形式为
)
x (f min x sub.to
0)x (C ≤0
)x (Ceq =b
x A ≤⋅beq
x Aeq =⋅0
)w ,x (K 11≤0
)w ,x (K 22≤…
)w ,x (K n n ≤其中：x 、b 、beq 、lb 、ub 都是向量；A 、Aeq 是矩阵；C(x)、Ceq(x)、)w ,x (K i i 是返回向量的函数，f(x)为目标函数；f(x)、C(x)、Ceq(x)是非线性函数；)w ,x (K i i 为半无限约束，n 21w ,,w ,w 通常是长度为2的向量。

在MTALAB5.x 中，使用函数seminf 解决这类问题。

187
函数fseminf
格式
x =fseminf(fun,x0,ntheta,seminfcon)
x =fseminf(fun,x0,ntheta,seminfcon,A,b)x =fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq)x =fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq,lb,ub)x =fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq,lb,ub,options)[x,fval]=fseminf(…)[x,fval,exitflag]=fseminf(…)[x,fval,exitflag,output]=fseminf(…)[x,fval,exitflag,output,lambda]=fseminf(…)
参数说明：x0为初始估计值；
fun 为目标函数，其定义方式与前面相同；A 、b 由线性不等式约束b x A ≤⋅确定，没有，则A=[]，b=[]；Aeq 、beq 由线性等式约束beq x Aeq =⋅确定，没有，则Aeq=[]，
beq=[]；
Lb 、ub 由变量x 的范围b u x b l ≤≤确定；
options 为优化参数；ntheta 为半无限约束的个数；seminfcon 用来确定非线性约束向量C 和Ceq 以及半无限约束的向
量K 1，K 2，…，K n ，通过指定函数柄来使用，如：x =fseminf(@myfun,x0,ntheta,@myinfcon)
先建立非线性约束和半无限约束函数文件，并保存为
myinfcon.m ：function [C,Ceq,K1,K2,…,Kntheta,S]=myinfcon(x,S)%S 为向量w 的采样值%初始化样本间距if isnan(S(1,1))，S =…%S 有ntheta 行2列end w1=…%计算样本集w2=…%计算样本集
…
wntheta =…%计算样本集K1=…%在x 和w 处的第1个半无限约束值K2=…%在x 和w 处的第2个半无限约束值
…
Kntheta =…%在x 和w 处的第ntheta 个半无限约束值
188
C =…%在x 处计算非线性不等式约束值Ceq =…%在x 处计算非线性等式约束值如果没有约束，则相应的值取为“[]”，如Ceq=[]fval 为在x 处的目标函数最小值；exitflag 为终止迭代的条件；output 为输出的优化信息；lambda 为解x 的Lagrange 乘子。

例5-10求下面一维情形的最优化问题
2
32221x )5.0x ()5.0x ()5.0x ()x (f min -+-+-=sub.to
1x )x w sin()50w (1000
1)x w cos()x w sin()w ,x (K 33121211111≤----=1x )x w sin()50w (1000
1)x w cos()x w sin()w ,x (K 33222122222≤----=100
w 11≤≤100
w 12≤≤解：将约束方程化为标准形式：
01x )x w sin()50w (1000
1)x w cos()x w sin()w ,x (K 33121211111≤-----=01x )x w sin()50w (1000
1)x w cos()x w sin()w ,x (K 33222122222≤-----=先建立非线性约束和半无限约束函数文件，并保存为mycon.m ：function [C,Ceq,K1,K2,S]=mycon(X,S)%初始化样本间距：if isnan(S(1,1)),S =[0.20;0.20];end %产生样本集：w1=1:S(1,1):100;w2=1:S(2,1):100;%计算半无限约束：K1=sin(w1*X(1)).*cos(w1*X(2))-1/1000*(w1-50).^2-sin(w1*X(3))-X(3)-1;K2=sin(w2*X(2)).*cos(w2*X(1))-1/1000*(w2-50).^2-sin(w2*X(3))-X(3)-1;%无非线性约束：C =[];Ceq=[];%绘制半无限约束图形plot(w1,K1,'-',w2,K2,':'),title('Semi-infinite constraints')
然后在MATLAB 命令窗口或编辑器中建立M 文件：fun ='sum((x-0.5).^2)';x0=[0.5;0.2;0.3];%Starting guess [x,fval]=fseminf(fun,x0,2,@mycon)
189
结果为：x =0.66730.30130.4023fval =0.0770>>[C,Ceq,K1,K2]=mycon (x,NaN);%利用初始样本间距
>>max(K1)ans =-0.0017>>max(K2)ans =-0.0845图5-1
例5-11求下面二维情形的最优化问题
2
32221x )2.0x ()2.0x ()2.0x ()x (f min -+-+-=sub.to
+----=3312122111x )x w sin()50w (1000
1)x w cos()x w sin()w ,x (K 5.1x )x w sin()50w (1000
1)x w cos()x w sin(332221122≤----100
w 11≤≤100
w 12≤≤初始点为x0=[0.25,0.25,0.25]。

解：先建立非线性和半无限约束函数文件，并保存为mycon.m ：function [C,Ceq,K1,S]=mycon(X,S)%初始化样本间距：if isnan(s(1,1)),s =[22];end %设置样本集w1x =1:s(1,1):100;w1y =1:s(1,2):100;[wx,wy]=meshgrid(w1x,w1y);%计算半无限约束函数值K1=sin(wx*X(1)).*cos(wx*X(2))-1/1000*(wx-50).^2
190
-sin(wx*X(3))-X(3)+…sin(wy*X(2)).*cos(wx*X(1))-1/1000*(wy-50).^2-sin(wy*X(3))-X(3)-1.5;%无非线性约束C =[];Ceq=[];%作约束曲面图形m =surf(wx,wy,K1,'edgecolor','none','facecolor','interp');camlight headlight title('Semi-infinite constraint')drawnow
然后在MATLAB 命令窗口下键入命令：>>fun ='sum((x-0.2).^2)';>>x0=[0.25,0.25,0.25];>>[x,fval]=fseminf(fun,x0,1,@mycon)
结果为（如图）x =0.29260.18740.2202fval =0.0091>>[c,ceq,K1]=mycon(x,[0.5,0.5]);%样本间距为0.5>>max(max(K1))ans =-0.0027
5.5极小化极大（Minmax ）问题
极小化极大问题的标准形式为
)}
x (F {max min i }F {x i sub.to 0
)x (C ≤0
)x (Ceq =b
x A ≤⋅beq
x Aeq =⋅ub
x lb ≤≤其中：x 、b 、beq 、lb 、ub 是向量，A 、Aeq 为矩阵，C(x)、Ceq(x)和F(x)是返回向量的函数，F(x)、C(x)、Ceq(x)可以是非线性函数。

在MATLAB5.x 中，它的求解由函数minmax 实现。

函数fminimax
格式x =fminimax(fun,x0)
x =fminimax(fun,x0,A,b)x =fminimax(fun,x0,A,b,Aeq,beq)x =
fminimax(fun,x0,A,b,Aeq,beq,lb,ub)x =fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
图5-2
191
x =fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)[x,fval,maxfval]=fminimax(…)[x,fval,maxfval,exitflag]=fminimax(…)[x,fval,maxfval,exitflag,output]=fminimax(…)[x,fval,maxfval,exitflag,output,lambda]=fminimax(…)
参数说明：fun 为目标函数；
x0为初始值；A 、b 满足线性不等约束b x A ≤⋅，若没有不等约束，则取A=[]，
b=[]；Aeq 、beq 满足等式约束beq x Aeq =⋅，若没有，则取Aeq=[]，beq=[]；
lb 、ub 满足ub x lb ≤≤，若没有界，可设lb=[]，ub=[]；
nonlcon 的作用是通过接受的向量x 来计算非线性不等约束0)x (C ≤和
等式约束0)x (Ceq =分别在x 处的值C 和Ceq ，通过指定函数柄来
使用，如：
>>x =fminimax(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非线性约束函数，并保存为mycon.m ：function [C,Ceq]=
mycon(x)
C =…
%计算x 处的非线性不等约束0)x (C ≤的函数值。

Ceq =…%计算x 处的非线性等式约束0)x (Ceq =的函数值。

options 为指定的优化参数；fval 为最优点处的目标函数值；maxfval 为目标函数在x 处的最大值；exitflag 为终止迭代的条件；lambda 是Lagrange 乘子，它体现哪一个约束有效。

output 输出优化信息。

例5-12求下列函数最大值的最小化问题
]
)x (f ,)x (f ,)x (f ,)x (f ,)x (f [54321其中：304
x 40x 48x x 2)x (f 2122211+--+=22
222x 3x )x (f --=18
x 3x )x (f 213-+=2
14x x )x (f --=8
x x )x (f 215-+=解：先建立目标函数文件，并保存为myfun.m ：function f =myfun(x)
f(1)=2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304;f(2)=-x(1)^2-3*x(2)^2;
192
f(3)=x(1)+3*x(2)-18;f(4)=-x(1)-x(2);f(5)=x(1)+x(2)-8;
然后，在命令窗口键入命令：x0=[0.1;0.1];%初始值[x,fval]=fminimax(@myfun,x0)
结果为：x = 4.00004.0000fval =0.0000-64.0000-2.0000-8.0000-0.0000
例5-13求上述问题的绝对值的最大值最小化问题。

目标函数为：]
|)x (f | ,|)x (f | ,|)x (f | ,|)x (f | ,|)x (f |[54321解：先建立目标函数文件（与上例相同）
然后，在命令窗口或编辑器中建立M 文件：>>x0=[0.1;0.1];%初始点>>options =optimset('MinAbsMax',5);%指定绝对值的最小化>>[x,fval]=fminimax(@myfun,x0,[],[],[],[],[],[],[],options)
则结果为x = 4.92562.0796fval =37.2356-37.2356-6.8357-7.0052-0.9948
5.6多目标规划问题
多目标规划是指在一组约束下，对多个不同目标函数进行优化。

它的一般形式为
]
)x (f ,),x (f ),x (f [min m 21 sub.to p
,,2,1j 0)x (g j =≤其中：)x ,,x ,x (x n 21 =。

在同一约束下，当目标函数处于冲突状态时，不存在最优解x 使所有目标函数同时达到最优。

此时，我们使用有效解，即如果不存在S x ∈，使得)x (f )x (f *i i ≥，i=1,2,…m,则称x*为有效解。

在MATLAB 中，多目标问题的标准形式为
γγ
,x imize min sub.to goal
weight )x (F ≤γ⋅-0
)x (C ≤
193
)x (Ceq =b
x A ≤⋅beq
x Aeq =⋅ub
x lb ≤≤其中：x 、b 、beq 、lb 、ub 是向量；A 、Aeq 为矩阵；C(x)、Ceq(x)和F(x)是返回向量的函数；F(x)、C(x)、Ceq(x)可以是非线性函数；weight 为权值系数向量，用于控制对应的目标函数与用户定义的目标函数值的接近程度；goal 为用户设计的与目标函数相应的目标函数值向量；γ为一个松弛因子标量；F(x)为多目标规划中的目标函数向量。

在MATLAB5.x 中，它的最优解由attgoal 函数实现。

函数fgoalattain
格式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 为多目标函数的文件名字符串，其定义方式与前面fun 的定义方
式相同；
goal 为用户设计的目标函数值向量；weight 为权值系数向量，用于控制目标函数与用户自定义目标值的接
近程度；
A 、b 满足线性不等式约束b x A ≤⋅，没有时取A=[]，b=[]；
Aeq 、beq 满足线性等式约束beq x Aeq =⋅，没有时取Aeq=[]，
beq=[]；lb 、ub 为变量的下界和上界：ub x lb ≤≤；
nonlcon 的作用是通过接受的向量x 来计算非线性不等约束0)x (C ≤和
等式约束0)x (Ceq =分别在x 处的值C 和Ceq ，通过指定函数柄
来使用。

194
如：
>>x =fgoalattain(@myfun,x0,goal,weight,A,b,Aeq,beq,lb,ub,@mycon)，先
建立非线性约束函数，并保存为mycon.m ：function [C,Ceq]=
mycon(x)
C =…
%计算x 处的非线性不等式约束0)x (C ≤的函数值。

Ceq =…
%计算x 处的非线性等式约束0)x (Ceq =的函数值。

options 为指定的优化参数；fval 为多目标函数在x 处的值；attainfactor 为解x 处的目标规划因子；exitflag 为终止迭代的条件；output 为输出的优化信息；lambda 为解x 处的Lagrange 乘子
例5-14控制系统输出反馈器设计。

设如下线性系统
Bu Ax x
+= Cx
y =其中：⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡---=210102000
5.o A ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=102201B ⎥⎦⎤⎢⎣⎡=100001C 要求设计输出反馈控制器K ，使闭环系统
u B x )BKC A (x
++= x C y =在复平面实轴上点[-5，-3，-1]的左侧有极点，并要求)2,1j ,i (4K 4j i =≤≤-解：上述问题就是要求解矩阵K ，使矩阵（A+BKC ）的极点为[-5，-3，-1]，这是一个多目标规划问题。

先建立目标函数文件，保存为eigfun.m ：
function F =eigfun(K,A,B,C)F =sort(eig(A+B*K*C));%估计目标函数值
然后，输入参数并调用优化程序：A =[-0.500;0-210;01-2];B =[10;-22;01];C =[100;001];K0=[-1-1;-1-1];%初始化控制器矩阵goal =[-5-3-1];%为闭合环路的特征值（极点）设置目标值向量weight =abs(goal)%设置权值向量
lb=-4*ones(size(K0));%设置控制器的下界
ub=4*ones(size(K0));%设置控制器的上界
options=optimset('Display','iter');%设置显示参数：显示每次迭代的输出
[K,fval,attainfactor]= fgoalattain(@eigfun,K0,goal,weight,[],[],[],[],lb,ub,[],options,A,B,C)
结果为：
weight=
531
Attainment Directional Iter F-count factor Step-size derivative Procedure
16 1.8851 1.03
213 1.0611-0.679
3200.42111-0.523 Hessian modified
427-0.063521-0.053 Hessian modified twice
534-0.15711-0.133
641-0.34891-0.00768 Hessian modified
748-0.36431-4.25e-005 Hessian modified
855-0.36451-0.00303 Hessian modified twice
962-0.36741-0.0213 Hessian modified
1069-0.380610.00266
1176-0.38621-2.73e-005 Hessian modified twice
1283-0.38631-1.22e-013 Hessian modified twice
Optimization terminated successfully:
Search direction less than2*options.TolX and maximum constraint violation is less than options.TolCon
Active Constraints:
1
2
4
9
10
K=
-4.0000-0.2564
-4.0000-4.0000
fval=
-6.9313
-4.1588
-1.4099
attainfactor=
-0.3863
195
196 5.7
最小二乘最优问题
5.7.1约束线性最小二乘有约束线性最小二乘的标准形式为
2
2
x d x C 21min -sub.to b x A ≤⋅beq
x Aeq =⋅b
u x b l ≤≤其中：C 、A 、Aeq 为矩阵；d 、b 、beq 、lb 、ub 、x 是向量。

在MATLAB5.x 中，约束线性最小二乘用函数conls 求解。

函数lsqlin
格式x =lsqlin(C,d,A,b)%求在约束条件b x A ≤⋅下，方程Cx =d 的最小二乘解x 。

x =lsqlin(C,d,A,b,Aeq,beq)%Aeq 、beq 满足等式约束beq x Aeq =⋅，
若没有不等式约束，则设A=[]，b=[]。

x =lsqlin(C,d,A,b,Aeq,beq,lb,ub)%lb 、ub 满足b u x b l ≤≤，若没有等
式约束，则Aeq=[]，beq=[]。

x =lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)%x0为初始解向量，若x 没有界，
则lb=[]，ub=[]。

x =lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options)%options 为指定优化参
数
[x,resnorm]=lsqlin(…)%resnorm=norm(C*x-d)^2，即2-范数。

[x,resnorm,residual]=lsqlin(…)%residual=C*x-d ，即残差。

[x,resnorm,residual,exitflag]=lsqlin(…)%exitflag 为终止迭代的条
件
[x,resnorm,residual,exitflag,output]=lsqlin(…)
%output 表示输出优
化信息
[x,resnorm,residual,exitflag,output,lambda]=lsqlin(…)%lambda 为解x 的Lagrange 乘子
例5-15求解下面系统的最小二乘解
系统：Cx=d
约束：b x A ≤⋅；b
u x b l ≤≤先输入系统系数和x 的上下界：C =[0.95010.76200.6153
0.4057;…
197
0.23110.45640.79190.9354;…0.60680.01850.92180.9169;…0.48590.82140.73820.4102;…0.89120.44470.17620.8936];d =[0.0578;0.3528;0.8131;0.0098;0.1388];A =[0.20270.27210.74670.4659;…0.19870.19880.44500.4186;…0.60370.01520.93180.8462];b =[0.5251;0.2026;0.6721];lb =-0.1*ones(4,1);ub =2*ones(4,1);
然后调用最小二乘命令：[x,resnorm,residual,exitflag,output,lambda]=lsqlin(C,d,A,b,[],[],lb,ub);结果为：x =-0.1000-0.10000.21520.3502resnorm =0.1672residual =0.04550.0764-0.35620.16200.0784exitflag =1%说明解x 是收敛的output =iterations:4algorithm:'medium-scale:active-set'firstorderopt:[]cgiterations:[]lambda =lower:[4x1double]upper:[4x1double]eqlin:[0x1double]ineqlin:[3x1double]
通过lambda.ineqlin 可查看非线性不等式约束是否有效。

5.7.2非线性数据（曲线）拟合
非线性曲线拟合是已知输入向量xdata 和输出向量ydata ，并且知道输入与输出的函数关系为ydata=F(x,xdata)，但不知道系数向量x 。

今进行曲线拟合，求x 使得下式成立：
∑-=-i
2i i 2
2x )ydata )xdata ,x (F (21ydata )xdata ,x (F 21min 在MATLAB5.x 中，使用函数curvefit 解决这类问题。

函数lsqcurvefit
格式x =lsqcurvefit(fun,x0,xdata,ydata)
198
x =lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x =lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)[x,resnorm]=lsqcurvefit(…)[x,resnorm,residual]=lsqcurvefit(…)[x,resnorm,residual,exitflag]=lsqcurvefit(…)[x,resnorm,residual,exitflag,output]=lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(…)
参数说明：
x0为初始解向量；xdata ，ydata 为满足关系ydata=F(x,xdata)的数据；lb 、ub 为解向量的下界和上界b u x b l ≤≤，若没有指定界，则lb=[]，
ub=[]；options 为指定的优化参数；fun 为拟合函数，其定义方式为：x =lsqcurvefit(@myfun,x0,xdata,ydata)，
其中myfun 已定义为
function F =myfun(x,xdata)F =…%计算x 处拟合函数值fun 的用法与前面相
同；
resnorm=sum ((fun(x,xdata)-ydata).^2)，即在x 处残差的平方和；residual=fun(x,xdata)-ydata ，即在x 处的残差；exitflag 为终止迭代的条件；output 为输出的优化信息；lambda 为解x 处的Lagrange 乘子；jacobian 为解x 处拟合函数fun 的jacobian 矩阵。

例5-16求解如下最小二乘非线性拟合问题
已知输入向量xdata 和输出向量ydata ，且长度都是n ，拟合函数为
3
2)i (xdata )3(x ))i (xdata sin()2(x )i (xdata )1(x )i (ydata ⋅+⋅+⋅=即目标函数为∑=-n
1i 2i i x )ydata )xdata ,x (F (21min 其中：3
2xdata )3(x )xdata sin()2(x xdata )1(x )xdata ,x (F ⋅+⋅+⋅=初始解向量为x0=[0.3,0.4,0.1]。

解：先建立拟合函数文件，并保存为myfun.m function F =myfun(x,xdata)F =x(1)*xdata.^2+x(2)*sin(xdata)+x(3)*xdata.^3;
然后给出数据xdata 和ydata >>xdata =[3.67.79.34.18.62.81.37.910.05.4];>>ydata =[16.5150.6263.124.7208.59.92.7163.9325.054.3];>>x0=[10,10,10];%初始估计值
199
>>[x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata)
结果为：Optimization terminated successfully:Relative function value changing by less than OPTIONS.TolFun x =0.22690.33850.3021resnorm =6.2950
5.7.3非线性最小二乘
非线性最小二乘（非线性数据拟合）的标准形式为
L
)x (f )x (f )x (f )x (f min 2m 2221x ++++= 其中：L 为常数
在MATLAB5.x 中，用函数leastsq 解决这类问题，在6.0版中使用函数lsqnonlin 。

设⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=)x (f )x (f )x (f )x (F m 21 则目标函数可表达为∑=i
2i 2
2x )x (f 21)x (F 21min 其中：x 为向量，F(x)为函数向量。

函数lsqnonlin
格式x =lsqnonlin(fun,x0)
%x0为初始解向量；fun 为)x (f i ，i=1,2,…,m ，fun 返回向量值F ，而不是平方和值，平方和隐含在算法中，fun
的定义与前面相同。

x =lsqnonlin(fun,x0,lb,ub)%lb 、ub 定义x 的下界和上界：
b u x b l ≤≤。

x =lsqnonlin(fun,x0,lb,ub,options)%options 为指定优化参数，若x 没
有界，则lb=[]，ub=[]。

[x,resnorm]=lsqnonlin(…)%resnorm=sum(fun(x).^2)，即解x 处目
标函数值。

[x,resnorm,residual]=lsqnonlin(…)%residual=fun(x)，即解x 处fun
的值。

[x,resnorm,residual,exitflag]=lsqnonlin(…)%exitflag 为终止迭代
条件。

[x,resnorm,residual,exitflag,output]=lsqnonlin(…)%output 输出优化
信息。

200
[x,resnorm,residual,exitflag,output,lambda]=lsqnonlin(…)%lambda 为Lagrage 乘子。

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(…)%fun 在解x 处的Jacobian 矩阵。

例5-17求下面非线性最小二乘问题∑=--+10
1k 2x k x k )e e k 22(21初始解向量为
x0=[0.3,0.4]。

解：先建立函数文件，并保存为myfun.m ，由于lsqnonlin 中的fun 为向量形式而不是平方和形式，因此，myfun 函数应由)x (f i 建立：
2
1x k x k k e e k 22)x (f --+=k=1,2,…,10function F =myfun(x)k =1:10;F =2+2*k-exp(k*x(1))-exp(k*x(2));
然后调用优化程序：x0=[0.30.4]；[x,resnorm]=lsqnonlin(@myfun,x0)
结果为：Optimization terminated successfully:Norm of the current step is less than OPTIONS.TolX x =0.25780.2578resnorm =%求目标函数值124.3622
5.7.4非负线性最小二乘非负线性最小二乘的标准形式为：
2
2
x d x C 21min -sub.to 0
x ≥其中：矩阵C 和向量d 为目标函数的系数，向量x 为非负独立变量。

在MATLAB5.x 中，用函数nnls 求解这类问题，在6.0版中则用函数lsqnonneg 。

函数lsqnonneg
格式x =lsqnonneg(C,d)%C 为实矩阵，d 为实向量
x =lsqnonneg(C,d,x0)%x0为初始值且大于0x =lsqnonneg(C,d,x0,options)%options 为指定优化参数[x,resnorm]=lsqnonneg(…)%resnorm=norm (C*x-d)^2[x,resnorm,residual]=lsqnonneg(…)%residual=C*x-d [x,resnorm,residual,exitflag]=lsqnonneg(…)[x,resnorm,residual,exitflag,output]=lsqnonneg(…)[x,resnorm,residual,exitflag,output,lambda]=lsqnonneg(…)
201
例5-18一个最小二乘问题的无约束与非负约束解法的比较。

先输入数据：>>C =[0.03720.2869;0.68610.7071;0.62330.6245;0.63440.6170];
>>d =[0.8587;0.1781;0.0747;0.8405];>>[C\d,lsqnonneg(C,d)]ans =-2.562703.11080.6929注意：1。

当问题为无约束线性最小二乘问题时，使用MATLAB 下的“\”运算即可以解决。

2．对于非负最小二乘问题，调用lsqnonneg(C,d)求解。

5.8
非线性方程(组)求解
5.8.1非线性方程的解
非线性方程的标准形式为f(x)=0
函数fzero
格式x =fzero (fun,x0)%用fun 定义表达式f(x)，x0为初始解。

x =fzero (fun,x0,options)[x,fval]=fzero(…)%fval=f(x)[x,fval,exitflag]=fzero(…)[x,fval,exitflag,output]=fzero(…)
说明该函数采用数值解求方程f(x)=0的根。

例5-19求05x 2x 3=--的根解：>>fun='x^3-2*x-5';>>z=fzero(fun,2)%初始估计值为2
结果为z = 2.0946
5.8.2非线性方程组的解
非线性方程组的标准形式为：F(x)=0
其中：x 为向量，F(x)为函数向量。

函数fsolve
格式x =fsolve(fun,x0)%用fun 定义向量函数，其定义方式为：先定义方
程函数function F =myfun (x)。

F =[表达式1；表达式2；…表达式m]%保存为myfun.m ，并用下面
方式调用：x =fsolve(@myfun,x0)，x0为初始估
计值。

x =fsolve(fun,x0,options)
202
[x,fval]=fsolve(…)%fval=F(x)，即函数值向量[x,fval,exitflag]=fsolve(…)[x,fval,exitflag,output]=fsolve(…)[x,fval,exitflag,output,jacobian]=fsolve(…)%jacobian 为解x 处的Jacobian 阵。

其余参数与前面参数相似。

例5-20求下列系统的根
2
1
x 21x 21e x 2x e x x 2--=+-=-解：化为标准形式
e x 2x 0
e x x 221x 21x 21=-+-=----设初值点为x0=[-5-5]。

先建立方程函数文件，并保存为myfun.m ：function F =myfun(x)F =[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];
然后调用优化程序x0=[-5;-5];%初始点options=optimset('Display','iter');%显示输出信息[x,fval]=fsolve(@myfun,x0,options)
结果为Norm of First-order Iteration Func-count f(x)step optimality CG-iterations 1447071.21 2.29e+0040276527.47 1.45207 3.09e+0031310918.372 1.491864181413127.74 1.5532657.3151614.9153 1.575918.2616190.779051 1.27662 1.1417220.003724530.4846580.068318259.21617e-0080.03855520.0003361928 5.66133e-0170.0001937078.34e-0091Optimization terminated successfully:Relative function value changing by less than OPTIONS.TolFun
203x =0.56710.5671fval =1.0e-008*-0.5320-0.5320
例5-21求矩阵x 使其满足方程⎥⎦⎤⎢⎣⎡=**4321x x x ，并设初始解向量为x=[1,1;1,1]。

解：先编写M 文件：
function F =myfun(x)F =x*x*x-[1,2;3,4];然后调用优化程序求解：>>x0=ones(2,2);%初始解向量
>>options =optimset('Display','off');%不显示优化信息>>[x,Fval,exitflag]=fsolve(@myfun,x0,options)则结果为x =-0.12910.8602
1.2903 1.1612
Fval =
1.0e-003*
0.1541-0.11630.0109-0.0243
exitflag =
1。