染雾matlab线性规划
MATLAB 优化问题
1 线性规划问题
线性规划问题是目标函数和约束条件均为线性函数的问题，MATLAB6.0解决的线性规划问题的标准形式为：
min f(x)x?R
sub.to：A?x?b
Aeq?x?beq
lb?x?ub n
其中f、x、b、beq、lb、ub为向量，A、Aeq为矩阵。
其它形式的线性规划问题都可经过适当变换化为此标准形式。
在MATLAB6.0版中，线性规划问题（Linear Programming）已用函数linprog取代了函数 linprog
格式 x = linprog(f,A,b) %求min f ' *x sub.to A?x?b线性规划的最优解。
x = linprog(f,A,b,Aeq,beq) %等式约束Aeq?x?beq，若没有不等式约束
A?x?b，则A=[ ]，b=[ ]。
x = linprog(f,A,b,Aeq,beq,lb,ub) %指定x的范围lb?x?ub，若没有等式约束
Aeq?x?beq ，则Aeq=[ ]，beq=[ ]
x = linprog(f,A,b,Aeq,beq,lb,ub,x0) %设置初值x0
x = 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
例5-1 求下面的优化问题
min ?5x1?4x2?6x3
sub.to x1?x2?x3?20
3x1?2x2?4x3?42
3x1?2x2?30
.
0?x1,0?x2,0?x3
解：
>>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)
结果为：
x = %最优解
0.0000
15.0000
3.0000
fval = %最优值
-78.0000
exitflag = %收敛
1
output =
iterations: 6 %迭代次数
cgiterations: 0
algorithm: 'lipsol' %所使用规则
lambda =
ineqlin: [3x1 double]
eqlin: [0x1 double]
upper: [3x1 double]
lower: [3x1 double]
>> lambda.ineqlin
ans =
0.0000
1.5000
0.5000
>> lambda.lower
ans =
1.0000
0.0000
0.0000
表明：不等约束条件2和3以及第1个下界是有效的
2 非线性规划问题
2.1 有约束的一元函数的最小值
单变量函数求最小值的标准形式为minf(x) sub.to x1?x?x2 x
在MATLAB中使用fmin函数求其最小值。
函数 fminbnd
格式 x = fminbnd(fun,x1,x2) %返回自变量x在区间x1?x?x2上函数fun取最小值
时x值，fun为目标函数的表达式字符串或MATLAB
自定义函数的函数柄。
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
例5-2 计算下面函数在区间(0，1)内的最小值。
x3?cosx?xlogx f(x)?e解：>> [x,fval,exitflag,output]=fminbnd('(x^3+cos(x)+x*log(x))/exp(x)',0,1)
x =
0.5223
fval =
0.3974
exitflag =
1
output =
iterations: 9
funcCount: 9
algorithm: 'golden section search, parabolic interpolation'
例5-3 在[0，5]上求下面函数的最小值f(x)?(x?3)3?1
解：先自定义函数：在MATLAB编辑器中建立M文件为：
function f = myfun(x)
f = (x-3).^2 - 1;
保存为myfun.m，然后在命令窗口键入命令：
>> x=fminbnd(@myfun,0,5)
则结果显示为：
x =
3
2.2 无约束多元函数最小值
多元函数最小值的标准形式为minf(x) x
其中：x为向量，如x?[x1,x2,?,xn]
在MATLAB中使用fmins求其最小值。
命令 利用函数fminsearch求无约束多元函数最小值
函数 fminsearch
格式 x = fminsearch(fun,x0) %x0为初始点，fun为目标函数的表达式字符串或
MATLAB自定义函数的函数柄。
x = fminsearch(fun,x0,options) % options查optimset
.
[x,fval] = fminsearch(?) %最优点的函数值
[x,fval,exitflag] = fminsearch(?) % exitflag与单变量情形一致
[x,fval,exitflag,output] = fminsearch(?) %output与单变量情形一致
注意：fminsearch采用了Nelder-Mead型简单搜寻法。
32例5-4 求y?2x1?4x1x32?10x1x2?x2的最小值点
解：>>X=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.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.0016 0.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效果较好。
2?2x1x2?x2例5-5 求f(x)?3x12的最小值。
>> fun='3*x(1)^2+2*x(1)*x(2)+x(2)^2';
>> x0=[1 1];
>> [x,fval,exitflag,output,grad,hessian]=fminunc(fun,x0)
结果为：
x =
1.0e-008 *
-0.7591 0.2665
fval =
1.3953e-016
exitflag =
1
output =
iterations: 3
funcCount: 16
stepsize: 1.2353
firstorderopt: 1.6772e-007
algorithm: 'medium-scale: Quasi-Newton line search'
grad =
1.0e-006 *
-0.1677
0.0114
hessian =
6.0000 2.0000
2.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=[1 1]；
>> x=fminunc(fun,x0)
x =
1.0e-008 *
-0.7591 0.2665
2.3 有约束的多元函数最小值
非线性有约束的多元函数的标准形式为：
minf(x) x
sub.to C(x)?0
Ceq(x)?0
A?x?b
Aeq?x?beq
lb?x?ub
其中：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)
[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满足线性不等式约束A?x?b，若没有不等式约束，则取A=[ ]，b=[ ]；
Aeq、beq满足等式约束Aeq?x?beq，若没有，则取Aeq=[ ]，beq=[ ]；
lb、ub满足lb?x?ub，若没有界，可设lb=[ ]，ub=[ ]；
nonlcon的作用是通过接受的向量x来计算非线性不等约束C(x)?0和等式
约束Ceq(x)?0分别在x处的估计C和Ceq，通过指定函数柄来使用，
如：>>x = fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非
线性约束函数，并保存为mycon.m：function [C,Ceq] = mycon(x)
C = ? % 计算x处的非线性不等约束C(x)?0的函数值。
Ceq = ? % 计算x处的非线性等式约束Ceq(x)?0的函数值。
lambda是Lagrange乘子，它体现哪一个约束有效。
output输出优化信息；
grad表示目标函数在x处的梯度；
hessian表示目标函数在x处的Hessiab值。
例5-6 求下面问题在初始点（0，1）处的最优解
2min x1?x22?x1x2?2x1?5x2
sub.to ?(x1?1)2?x2?0
2x1?3x2?6?0
解：约束条件的标准形式为
sub.to (x1?1)2?x2?0
?2x1?3x2?6
先在MATLAB编辑器中建立非线性约束函数文件：
function [c, ceq]=mycon (x)
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=[0 1];
>>A=[-2 3]; %线性不等式约束
>>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 =
3 4
fval =
-13
exitflag = %解收敛
1
output =
iterations: 2
funcCount: 9
stepsize: 1
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: [ ]
cgiterations: [ ]
lambda =
lower: [2x1 double] %x下界有效情况，通过lambda.lower可查看。
upper: [2x1 double] %x上界有效情况，为0表示约束无效。
eqlin: [0x1 double] %线性等式约束有效情况，不为0表示约束有效。
eqnonlin: [0x1 double] %非线性等式约束有效情况。
ineqlin: 2.5081e-008 %线性不等式约束有效情况。
ineqnonlin: 6.1938e-008 %非线性不等式约束有效情况。
grad = %目标函数在最小值点的梯度
1.0e-006 *
-0.1776
hessian = %目标函数在最小值点的Hessian值
1.0000 -0.0000
-0.0000 1.0000
例5-7 求下面问题在初始点x=(10, 10, 10)处的最优解。
Min f(x)??x1x2x3
Sub.to 0?x1?2x2?2x3?72
解：约束条件的标准形式为
sub.to ?x1?2x2?2x3?0 x1?2x2?2x3?72
>> 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)
结果为：
x =
24.0000 12.0000 12.0000
fval =
-3456
2.4 二次规划问题
二次规划问题（quadratic programming）的标准形式为：
minx?Hx?f?x sub.to A?x?b
?x?beq Aeq
lb?x?ub
其中，H、A、Aeq为矩阵，f、b、beq、lb、ub、x为向量
.
其它形式的二次规划问题都可转化为标准形式。
MATLAB中二次规划应用函数quadprog。
函数 quadprog
格式 x = quadprog(H,f,A,b) %其中H,f,A,b为标准形中的参数，x为目标函数的最小
值。
x = quadprog(H,f,A,b,Aeq,beq) %Aeq,beq满足等约束条件Aeq?x?beq。
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 求解下面二次规划问题
min
sub.to 2f(x)?x1?x22?x1x2?2x1?6x2 x1?x2?2
?x1?2x2?2
2x1?x2?3
0?x1,0?x2
解：f(x)?x?Hx?f?x ?1?1??x1???2?f?x?则H??，，???6??x? ?12???2???
在MATLAB中实现如下：
>>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,b,[ ],[ ],lb)
结果为：
x = %最优解
0.6667
1.3333
fval = %最优值
-8.2222
exitflag = %收敛
1
output =
iterations: 3
algorithm: 'medium-scale: active-set'
firstorderopt: [ ]
cgiterations: [ ]
lambda =
lower: [2x1 double]
upper: [2x1 double]
eqlin: [0x1 double]
ineqlin: [3x1 double]
>> lambda.ineqlin
ans =
3.1111
0.4444
>> lambda.lower
ans =
说明 第1、2个约束条件有效，其余无效。
例5-9 求二次规划的最优解
max f (x1, x2)=x1x2+3
sub.to x1+x2-2=0
解：化成标准形式：
?0?1??x1??x1?minf(x1x2)??x1x2?3?(x1x2)?????(0,0)???3 ??10??x2??x2?
sub.to x1+x2=2
在Matlab中实现如下：
>>H=[0,-1;-1,0];
>>f=[0;0];
>>Aeq=[1 1];
>>b=2;
>>[x,fval,exitflag,output,lambda] = quadprog(H,f,[ ],[ ],Aeq,b)
结果为：
x =
1.0000
1.0000
fval =
-1.0000
exitflag =
1
output =
firstorderopt: 0
iterations: 1
cgiterations: 1
algorithm: [1x58 char]
lambda =
eqlin: 1.0000
ineqlin: [ ]
lower: [ ]
upper: [ ]
.
3 “半无限”有约束的多元函数最优解
“半无限”有约束多元函数最优解问题的标准形式为
minxf(x)
sub.to C(x)?0
Ceq(x)?0
A?x?b
Aeq?x?beq
K1(x,w1)?0
K2(x,w2)?0
?
Kn(x,wn)?0
其中：x、b、beq、lb、ub都是向量；A、Aeq是矩阵；C(x)、Ceq(x)、Ki(x,wi)是返回向量的函数，f(x)为目标函数；f(x)、C(x)、Ceq(x)是非线性函数；Ki(x,wi)为半无限约束，w1,w2,?,wn通常是长度为2的向量。
在MTALAB5.x中，使用函数seminf解决这类问题。
函数 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由线性不等式约束A?x?b确定，没有，则A=[ ]，b=[ ]；
Aeq、beq由线性等式约束Aeq?x?beq确定，没有，则Aeq=[ ]，beq=[ ]； Lb、ub由变量x的范围lb?x?ub确定；
options为优化参数；
ntheta为半无限约束的个数；
seminfcon用来确定非线性约束向量C和Ceq以及半无限约束的向量K1，
K2，?，Kn，通过指定函数柄来使用，如：
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个半无限约束值
C = ? % 在x处计算非线性不等式约束值
Ceq = ? % 在x处计算非线性等式约束值
如果没有约束，则相应的值取为“[ ]”，如Ceq=[]
fval为在x处的目标函数最小值；
exitflag为终止迭代的条件；
output为输出的优化信息；
lambda为解x的Lagrange乘子。
例5-10 求下面一维情形的最优化问题
minxf(x)?(x1?0.5)2?(x2?0.5)2?(x3?0.5)2
sub.to
K1(x,w1)?sin(w1x1)cos(w1x2)?(w1?50)2?sin(w1x3)?x3?1 K2(x,w2)?sin(w2x2)cos(w2x1)?(w2?50)2?sin(w2x3)?x3?1 1?w1?100
1?w2?100
解：将约束方程化为标准形式：
K1(x,w1)?sin(w1x1)cos(w1x2)?(w1?50)2?sin(w1x3)?x3?1?0 K2(x,w2)?sin(w2x2)cos(w2x1)?(w2?50)2?sin(w2x3)?x3?1?0 先建立非线性约束和半无限约束函数文件，并保存为mycon.m：
function [C,Ceq,K1,K2,S] = mycon(X,S)
% 初始化样本间距：
.
if isnan(S(1,1)),
S = [0.2 0; 0.2 0];
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)
结果为：
x =
0.6673
0.3013
0.4023
fval =
0.0770
>>[C,Ceq,K1,K2] = mycon (x,NaN); % 利用初始样本间距
>>max(K1)
ans =
-0.0017
>>max(K2)
ans =
-0.0845
图5-1
例5-11 求下面二维情形的最优化问题
minxf(x)?(x1?0.2)2?(x2?0.2)2?(x3?0.2)2
sub.to
K1(x,w)?sin(w1x1)cos(w2x2)?(w1?50)2?sin(w1x3)?x3??
sin(w2x2)cos(w1x1)?(w2?50)2?sin(w2x3)?x3?1.5 1?w1?100
1?w2?100
初始点为x0=[0.25, 0.25, 0.25]。
解：先建立非线性和半无限约束函数文件，并保存为mycon.m：
function [C,Ceq,K1,S] = mycon(X,S)
% 初始化样本间距：
if isnan(s(1,1)),
s = [2 2];
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 -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.2926 0.1874 0.2202
fval =
0.0091
>>[c,ceq,K1] = mycon(x,[0.5,0.5]); % 样本间距为0.5
>>max(max(K1))
ans =
-0.0027 图5-2
4 极小化极大（Minmax）问题
极小化极大问题的标准形式为
minmaxx{Fi}{Fi(x)}
sub.to C(x)?0
Ceq(x)?0
A?x?b
.
Aeq?x?beq
lb?x?ub
其中：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)
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满足线性不等约束A?x?b，若没有不等约束，则取A=[ ]，b=[ ]；
Aeq、beq满足等式约束Aeq?x?beq，若没有，则取Aeq=[ ]，beq=[ ]；
lb、ub满足lb?x?ub，若没有界，可设lb=[ ]，ub=[ ]；
nonlcon的作用是通过接受的向量x来计算非线性不等约束C(x)?0和等式约束
Ceq(x)?0分别在x处的值C和Ceq，通过指定函数柄来使用，如：>>x =
fminimax(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)，先建立非线性约束函数，
并保存为mycon.m：function [C,Ceq] = mycon(x)
C = ? % 计算x处的非线性不等约束C(x)?0的函数值。
Ceq = ? % 计算x处的非线性等式约束Ceq(x)?0的函数值。
options为指定的优化参数；
fval为最优点处的目标函数值；
maxfval为目标函数在x处的最大值；
exitflag为终止迭代的条件；
lambda是Lagrange乘子，它体现哪一个约束有效。
output输出优化信息。
例5-12 求下列函数最大值的最小化问题
[f1(x), f2(x), f3(x), f4(x), f5(x)]
2?x2其中：f1(x)?2x12?48x1?40x2?304
2f2(x)??x22?3x2
f3(x)?x1?3x2?18
f4(x)??x1?x2
f5(x)?x1?x2?8
解：先建立目标函数文件，并保存为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;
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.0000
4.0000
fval =
0.0000 -64.0000 -2.0000 -8.0000 -0.0000
例5-13 求上述问题的绝对值的最大值最小化问题。
目标函数为：[|f1(x)|, |f2(x)|, |f3(x)|, |f4(x)|, |f5(x)|]
解：先建立目标函数文件（与上例相同）
然后，在命令窗口或编辑器中建立M文件：
>>x0 = [0.1; 0.1]; % 初始点
>>options = optimset('MinAbsMax',5); % 指定绝对值的最小化
>>[x,fval] = fminimax(@myfun,x0,[ ],[ ],[ ],[ ],[ ],[ ],[ ],options)
则结果为
x =
4.9256
2.0796
fval =
37.2356 -37.2356 -6.8357 -7.0052 -0.9948
5 多目标规划问题
多目标规划是指在一组约束下，对多个不同目标函数进行优化。它的一般形式为
min[f1(x),f2(x),?,fm(x)]
j?1,2,?,p sub.to gj(x)?0
其中：x?(x1,x2,?,xn)。
在同一约束下，当目标函数处于冲突状态时，不存在最优解x使所有目标函数同时达到最优。此时，我们使用有效解，即如果不存在x?S，使得fi(x)?fi(x*)，i=1,2,?m, 则称x*为有效解。
在MATLAB中，多目标问题的标准形式为
.
minimize? x,?
sub.to F(x)?weight???goal
C(x)?0
Ceq(x)?0
A?x?b
Aeq?x?beq
lb?x?ub
其中：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满足线性不等式约束A?x?b，没有时取A=[ ]，b=[ ]；
Aeq、beq满足线性等式约束Aeq?x?beq，没有时取Aeq=[ ]，beq=[ ]；
lb、ub为变量的下界和上界：lb?x?ub；
nonlcon的作用是通过接受的向量x来计算非线性不等约束C(x)?0和等式约束
Ceq(x)?0分别在x处的值C和Ceq，通过指定函数柄来使用。
如：>>x = fgoalattain(@myfun,x0,goal,wei
/news/559EE48D172564E3.html ght,A,b,Aeq,beq,lb,ub,@mycon)，先建立非线性约束函数，并保存为mycon.m：function [C,Ceq] = mycon(x)
C = ? % 计算x处的非线性不等式约束C(x)?0的函数值。
Ceq = ? % 计算x处的非线性等式约束Ceq(x)?0的'函数值。
options为指定的优化参数；
fval为多目标函数在x处的值；
attainfactor为解x处的目标规划因子；
exitflag为终止迭代的条件；
output为输出的优化信息；
lambda为解x处的Lagrange乘子
例5-14 控制系统输出反馈器设计。
设如下线性系统
??Ax?Bu x
y?Cx
0???o.50?10??100??210? B???22? C??其中：A??0? ????001????1?2??0??01??
要求设计输出反馈控制器K，使闭环系统
??(A?BKC)x?Bu x
y?Cx
在复平面实轴上点[-5，-3，-1]的左侧有极点，并要求 ?4?Kij?4(i,j?1,2)
解：上述问题就是要求解矩阵K，使矩阵（A+BKC）的极点为[-5，-3，-1]，这是一个多目标规划问题。
先建立目标函数文件，保存为eigfun.m：
function F = eigfun(K,A,B,C)
F = sort(eig(A+B*K*C)); % 估计目标函数值
然后，输入参数并调用优化程序：
A = [-0.5 0 0; 0 -2 10; 0 1 -2];
B = [1 0; -2 2; 0 1];
C = [1 0 0; 0 0 1];
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 =
5 3 1
Attainment Directional
Iter F-count factor Step-size derivative Procedure
1 6 1.885 1 1.03
2 13 1.061 1 -0.679
3 20 0.4211 1 -0.523 Hessian modified
4 27 -0.06352 1 -0.053 Hessian modified twice 5 34 -0.1571 1 -0.133
.
6 41 -0.3489 1 -0.00768 Hessian modified
7 48 -0.3643 1 -4.25e-005 Hessian modified
8 55 -0.3645 1 -0.00303 Hessian modified twice 9 62 -0.3674 1 -0.0213 Hessian modified
10 69 -0.3806 1 0.00266
11 76 -0.3862 1 -2.73e-005 Hessian modified twice 12 83 -0.3863 1 -1.22e-013 Hessian modified twice Optimization terminated successfully:
Search direction less than 2*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
6 最小二乘最优问题
6.1 约束线性最小二乘
有约束线性最小二乘的标准形式为
minCx?dxsub.to A?x?b
Aeq?x?beq
lb?x?ub 22
其中：C、A、Aeq为矩阵；d、b、beq、lb、ub、x是向量。
在MATLAB5.x中，约束线性最小二乘用函数conls求解。
函数 lsqlin
格式 x = lsqlin(C,d,A,b) %求在约束条件A?x?b下，方程Cx = d的最小二乘解x。
x = lsqlin(C,d,A,b,Aeq,beq) %Aeq、beq满足等式约束Aeq?x?beq，若没有不
等式约束，则设A=[ ]，b=[ ]。
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub) %lb、ub满足lb?x?ub，若没有等式约束，
则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
约束：A?x?b；lb?x?ub
先输入系统系数和x的上下界：
C = [0.9501 0.7620 0.6153 0.4057;…
0.2311 0.4564 0.7919 0.9354;…
0.6068 0.0185 0.9218 0.9169;…
0.4859 0.8214 0.7382 0.4102;…
0.8912 0.4447 0.1762 0.8936];
d = [ 0.0578; 0.3528; 0.8131; 0.0098; 0.1388];
A =[ 0.2027 0.2721 0.7467 0.4659;…
0.1987 0.1988 0.4450 0.4186;…
0.6037 0.0152 0.9318 0.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.1000
0.2152
0.3502
resnorm =
0.1672
residual =
0.0455
0.0764
-0.3562
0.1620
0.0784
exitflag =
1 %说明解x是收敛的
output =
iterations: 4
algorithm: 'medium-scale: active-set'
firstorderopt: []
cgiterations: []
lambda =
lower: [4x1 double]
upper: [4x1 double]
.
eqlin: [0x1 double]
ineqlin: [3x1 double]
通过lambda.ineqlin可查看非线性不等式约束是否有效。
6.2 非线性数据（曲线）拟合
非线性曲线拟合是已知输入向量xdata和输出向量ydata，并且知道输入与输出的函数关系为ydata=F(x, xdata)，但不知道系数向量x。今进行曲线拟合，求x使得下式成立：
minF(x,xdata)?ydatax2
2??(F(x,xdatai)?ydatai)2 i
在MATLAB5.x中，使用函数curvefit解决这类问题。
函数 lsqcurvefit
格式 x = lsqcurvefit(fun,x0,xdata,ydata)
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为解向量的下界和上界lb?x?ub，若没有指定界，则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，拟合函数为
ydata(i)?x(1)?xdata(i)2?x(2)?sin(xdata(i))?x(3)?xdata(i)3
即目标函数为min?(F(x,xdatai)?ydatai)2 xi?1
其中：F(x,xdata)?x(1)?xdata2?x(2)?sin(xdata)?x(3)?xdata3
n
初始解向量为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.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4];
>>ydata = [16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3];
>>x0 = [10, 10, 10]; %初始估计值
>>[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata)
结果为：
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
x =
0.2269 0.3385 0.3021
resnorm =
6.2950
6.3 非线性最小二乘
非线性最小二乘（非线性数据拟合）的标准形式为
minf(x)?f1(x)2?f2(x)2???fm(x)2?L x
其中：L为常数
在MATLAB5.x中，用函数leastsq解决这类问题，在6.0版中使用函数lsqnonlin。
?f1(x)??f(x)?2? 设F(x)???????f(x)m??
则目标函数可表达为minF(x)x2
2??fi(x)2 i
其中：x为向量，F(x)为函数向量。
函数 lsqnonlin
格式 x = lsqnonlin(fun,x0) %x0为初始解向量；fun为fi(x)，i=1,2,?,m，fun返回向
量值F，而不是平方和值，平方和隐含在算法中，fun的定义与前面相同。
x = lsqnonlin(fun,x0,lb,ub) %lb、ub定义x的下界和上界：lb?x?ub。
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输出优化信息。
[x,resnorm,residual,exitflag,output,lambda] = lsqnonlin(?) %lambda为Lagrage
.
乘子。
[x,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(?) %fun在解x
处的Jacobian矩阵。
例5-17 求下面非线性最小二乘问题?(2?2k?ekx1?ekx2)2初始解向量为x0=[0.3,
k?110
0.4]。
解：先建立函数文件，并保存为myfun.m，由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由fi(x)建立：
fk(x)?2?2k?ekx1?ekx2 k=1,2,…,10
function F = myfun(x)
k = 1:10;
F = 2 + 2*k-exp(k*x(1))-exp(k*x(2));
然后调用优化程序：
x0 = [0.3 0.4]；
[x,resnorm] = lsqnonlin(@myfun,x0)
结果为：
Optimization terminated successfully:
Norm of the current step is less than OPTIONS.TolX
x =
0.2578 0.2578
resnorm = %求目标函数值
124.3622
6.4 非负线性最小二乘
非负线性最小二乘的标准形式为：
minCx?dxsub.to x?0 22
其中：矩阵C和向量d为目标函数的系数，向量x为非负独立变量。
在MATLAB5.x中，用函数nnls求解这类问题，在6.0版中则用函数lsqnonneg。 函数 lsqnonneg
格式 x = lsqnonneg(C,d) %C为实矩阵，d为实向量
x = lsqnonneg(C,d,x0) % x0为初始值且大于0
x = 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(?)
例5-18 一个最小二乘问题的无约束与非负约束解法的比较。
先输入数据：
>>C = [ 0.0372 0.2869; 0.6861 0.7071; 0.6233 0.6245; 0.6344 0.6170];
>>d = [0.8587; 0.1781; 0.0747; 0.8405];
>> [C\d, lsqnonneg(C,d)]
ans =
-2.5627 0
3.1108 0.6929
注意：1。当问题为无约束线性最小二乘问题时，使用MATLAB下的“\”运算即可以解决。2．对于非负最小二乘问题，调用lsqnonneg(C,d)求解。
7 非线性方程(组)求解
7.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 求x3?2x?5?0的根
解：>> fun='x^3-2*x-5';
>> z=fzero(fun,2) %初始估计值为2
结果为
z =
2.0946
7.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)
[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 求下列系统的根
2x1?x2?e?x1
?x1?2x2?e?x2
解：化为标准形式
2x1?x2?e?x1?0
?x1?2x2?e?x2?0
设初值点为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 1 4 47071.2 1 2.29e+004 0
2 7 6527.47 1.45207 3.09e+003 1
3 10 918.372 1.49186 418 1
4 13 127.74 1.55326 57.3 1
5 16 14.9153 1.57591 8.26 1
6 19 0.779051 1.27662 1.14 1
7 22 0.00372453 0.484658 0.0683 1
8 25 9.21617e-008 0.0385552 0.000336 1
9 28 5.66133e-017 0.000193707 8.34e-009 1
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
x =
0.5671
0.5671
fval =
1.0e-008 *
-0.5320
-0.5320
?12?例5-21 求矩阵x使其满足方程x?x?x???，并设初始解向量为x=[1, 1; 1, 1]。 34??
解：先编写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.1291 0.8602
1.2903 1.1612
Fval =
1.0e-003 *
0.1541 -0.1163
0.0109
exitflag =
1
. -0.0243
