MATLAB优化应用.docx
《MATLAB优化应用.docx》由会员分享,可在线阅读,更多相关《MATLAB优化应用.docx(22页珍藏版)》请在冰点文库上搜索。
MATLAB优化应用
MATLAB优化应用
§1线性规划模型
一、线性规划课题:
实例1:
生产计划问题
假设某厂计划生产甲、乙两种产品,现库存主要材料有A类3600公斤,B类2000公斤,C类3000公斤。
每件甲产品需用材料A类9公斤,B类4公斤,C类3公斤。
每件乙产品,需用材料A类4公斤,B类5公斤,C类10公斤。
甲单位产品的利润70元,乙单位产品的利润120元。
问如何安排生产,才能使该厂所获的利润最大。
建立数学模型:
设x1、x2分别为生产甲、乙产品的件数。
f为该厂所获总润。
maxf=70x1+120x2
s.t9x1+4x2≤3600
4x1+5x2≤2000
3x1+10x2≤3000
x1,x2≥0
实例2:
投资问题
某公司有一批资金用于4个工程项目的投资,其投资各项目时所得的净收益(投入资金锪百分比)如下表:
工程项目收益表
工程项目
A
B
C
D
收益(%)
15
10
8
12
由于某种原因,决定用于项目A的投资不大于其他各项投资之和而用于项目B和C的投资要大于项目D的投资。
试确定全文该公司收益最大的投资分配方案。
建立数学模型:
设x1、x2、x3、x4分别代表用于项目A、B、C、D的投资百分数。
maxf=0.15x1+0.1x2+0.08x3+0.12x4
s.tx1-x2-x3-x4≤0
x2+x3-x4≥0
x1+x2+x3+x4=1
xj≥0j=1,2,3,4
实例3:
运输问题
有A、B、C三个食品加工厂,负责供给甲、乙、丙、丁四个市场。
三个厂每天生产食品箱数上限如下表:
工厂
A
B
C
生产数
60
40
50
四个市场每天的需求量如下表:
市场
甲
乙
丙
丁
需求量
20
35
33
34
从各厂运到各市场的运输费(元/每箱)由下表给出:
市场
甲
乙
丙
丁
工
厂
A
2
1
3
2
B
1
3
2
1
C
3
4
1
1
求在基本满足供需平衡的约束条件下使总运输费用最小。
建立数学模型:
设aij为由工厂i运到市场j的费用,xij是由工厂i运到市场j的箱数。
bi是工厂i的产量,dj是市场j的需求量。
b=(604050)d=(20353334)
s.t
xij≥0
当我们用MATLAB软件作优化问题时,所有求maxf的问题化为求min(-f)来作。
约束gi(x)≥0,化为–gi≤0来作。
上述实例去掉实际背景,归结出规划问题:
目标函数和约束条件都是变量x的线性函数。
形如:
(1)minfTX
s.tAX≤b
AeqX=beq
lb≤X≤ub
其中X为n维未知向量,fT=[f1,f2,…fn]为目标函数系数向量,小于等于约束系数矩阵A为m×n矩阵,b为其右端m维列向量,Aeq为等式约束系数矩阵,beq为等式约束右端常数列向量。
lb,ub为自变量取值上界与下界约束的n维常数向量。
二.线性规划问题求最优解函数:
调用格式:
x=linprog(f,A,b)
x=linprog(f,A,b,Aeq,beq)
x=linprog(f,A,b,Aeq,beq,lb,ub)
x=linprog(f,A,b,Aeq,beq,lb,ub,x0)
x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)
[x,fval]=linprog(…)
[x,fval,exitflag]=linprog(…)
[x,fval,exitflag,output]=linprog(…)
[x,fval,exitflag,output,lambda]=linprog(…)
说明:
x=linprog(f,A,b)返回值x为最优解向量。
x=linprog(f,A,b,Aeq,beq)作有等式约束的问题。
若没有不等式约束,则令A=[]、b=[]。
x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)中lb,ub为变量x的下界和上界,x0为初值点,options为指定优化参数进行最小化。
Options的参数描述:
Display显示水平。
选择’off’不显示输出;选择’iter’显示每一步迭代过程的输出;选择’final’显示最终结果。
MaxFunEvals函数评价的最大允许次数
Maxiter最大允许迭代次数
TolXx处的终止容限
[x,fval]=linprog(…)左端fval返回解x处的目标函数值。
[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb,ub,x0)的输出部分:
exitflag描述函数计算的退出条件:
若为正值,表示目标函数收敛于解x处;若为负值,表示目标函数不收敛;若为零值,表示已经达到函数评价或迭代的最大次数。
output返回优化信息:
output.iterations表示迭代次数;output.algorithm表示所采用的算法;outprt.funcCount表示函数评价次数。
lambda返回x处的拉格朗日乘子。
它有以下属性:
lambda.lower-lambda的下界;
lambda.upper-lambda的上界;
lambda.ineqlin-lambda的线性不等式;
lambda.eqlin-lambda的线性等式。
三.举例
例1:
求解线性规划问题:
maxf=2x1+5x2
s.t
先将目标函数转化成最小值问题:
min(-f)=-2x1-5x2
程序:
f=[-2-5];
A=[10;01;12];
b=[4;3;8];
[x,fval]=linprog(f,A,b)
f=fval*(-1)
结果:
x=2
3
fval=-19.0000
maxf=19
例2:
minf=5x1-x2+2x3+3x4-8x5
s.t–2x1+x2-x3+x4-3x5≤6
2x1+x2-x3+4x4+x5≤7
0≤xj≤15j=1,2,3,4,5
程序:
f=[5-123-8];
A=[-21-11-3;21-141];
b=[6;7];
lb=[00000];
ub=[1515151515];
[x,fval]=linprog(f,A,b,[],[],lb,ub)
结果:
x=
0.0000
0.0000
8.0000
0.0000
15.0000
minf=
-104
例3:
求解线性规划问题:
minf=5x1+x2+2x3+3x4+x5
s.t–2x1+x2-x3+x4-3x5≤1
2x1+3x2-x3+2x4+x5≤-2
0≤xj≤1j=1,2,3,4,5
程序:
f=[51231];
A=[-21-11-3;23-121];
b=[1;-2];
lb=[00000];
ub=[11111];
[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb,ub)运行结果:
Exiting:
Oneormoreoftheresiduals,dualitygap,ortotalrelativeerror
hasgrown100000timesgreaterthanitsminimumvaluesofar:
theprimalappearstobeinfeasible(andthedualunbounded).
(Thedualresidual
x=0.0000
0.0000
1.1987
0.0000
0.0000
fval=
2.3975
exitflag=
-1
output=
iterations:
7
cgiterations:
0
algorithm:
'lipsol'
lambda=
ineqlin:
[2x1double]
eqlin:
[0x1double]
upper:
[5x1double]
lower:
[5x1double]
显示的信息表明该问题无可行解。
所给出的是对约束破坏最小的解。
例4:
求解实例1的生产计划问题
建立数学模型:
设x1、x2分别为生产甲、乙产品的件数。
f为该厂所获总润。
maxf=70x1+120x2
s.t9x1+4x2≤3600
4x1+5x2≤2000
3x1+10x2≤3000
x1,x2≥0
将其转换为标准形式:
minf=-70x1-120x2
s.t9x1+4x2≤3600
4x1+5x2≤2000
3x1+10x2≤3000
x1,x2≥0
程序:
f=[-70-120];
A=[94;45;310];
b=[3600;2000;3000];
lb=[00];
ub=[];
[x,fval,exitflag]=linprog(f,A,b,[],[],lb,ub)
maxf=-fval
结果:
x=
200.0000
240.0000
fval=
-4.2800e+004
exitflag=
1
maxf=
4.2800e+004
例5:
求解实例2
建立数学模型:
maxf=0.15x1+0.1x2+0.08x3+0.12x4
s.tx1-x2-x3-x4≤0
x2+x3-x4≥0
x1+x2+x3+x4=1
xj≥0j=1,2,3,4
将其转换为标准形式:
minz=-0.15x1-0.1x2-0.08x3-0.12x4
s.tx1-x2-x3-x4≤0
-x2-x3+x4≤0
x1+x2+x3+x4=1
xj≥0j=1,2,3,4
程序:
f=[-0.15;-0.1;-0.08;-0.12];
A=[1-1-1-1
0-1-11];
b=[0;0];
Aeq=[1111];
beq=[1];
lb=zeros(4,1);
[x,fval,exitflag]=linprog(f,A,b,Aeq,beq,lb)
f=-fval
结果:
x=
0.5000
0.2500
0.0000
0.2500
fval=
-0.1300
exitflag=
1
f=
0.1300
即4个项目的投资百分数分别为50%,25%,0,25%时可使该公司获得最大的收益,其最大收益可到达13%。
过程正常收敛。
例6:
求解实例3
建立数学模型:
设aij为由工厂i运到市场j的费用,xij是由工厂i运到市场j的箱数。
bi是工厂i的产量,dj是市场j的需求量。
b=(604050)Td=(20353334)T
s.t
xij≥0
程序:
A=[2132;1321;3411];
f=A(:
);
B=[100100100100
010010010010
001001001001];
D=[111000000000
000111000000
000000111000
000000000111];
b=[60;40;50];
d=[20;35;33;34];
lb=zeros(12,1);
[x,fval,exitflag]=linprog(f,B,b,D,d,lb)
结果:
x=
0.0000
20.0000
0.0000
35.0000
0.0000
0.0000
0.0000
0.0000
33.0000
0.0000
18.4682
15.5318
fval=
122.0000
exitflag=
1
即运输方案为:
甲市场的货由B厂送20箱;乙市场的货由A厂送35箱;丙商场的货由C厂送33箱;丁市场的货由B厂送18箱,再由C厂送16箱。
最低总运费为:
122元。
§2非线性规划模型
一.非线性规划课题
实例1表面积为36平方米的最大长方体体积。
建立数学模型:
程序:
首先建立目标函数文件ff8.m文件:
functionf=ff8(x)
f=exp(x
(1))*(6*x
(1)^2+3*x
(2)^2+2*x
(1)*x
(2)+4*x
(2)+1);
再建立非线性的约束条件文件:
ff8g.m
function[c,g]=ff8g(x)
c
(1)=x
(1)*x
(2)-x
(1)-x
(2)+1;
c
(2)=-2*x
(1)*x
(2)-5;
g=[];
然后在工作空间键入程序:
x0=[1,1];
nonlcon=@ff8g
[x,fval]=fmincon(@ff8,x0,[],[],[],[],[],[],nonlcon)
结果:
x=
-2.50001.0000
fval=
3.3244
exitflag=
1
当有等式约束时,要放在矩阵g的位置,如上例中加等式约束:
x
(1)+2*x
(1)=0
程序:
首先建立fun1.m文件:
function[c,g]=ff8g1(x)
c
(1)=x
(1)*x
(2)-x
(1)-x
(2)+1;
c
(2)=-2*x
(1)*x
(2)-5;
g
(1)=x
(1)+2*x
(2);
然后在工作空间键入程序:
x0=[-1,1];
nonlcon=@ff8g1;
[x,fval,exitflag]=fmincon(@ff8,x0,[],[],[],[],[],[],nonlcon)
结果:
x=
-2.23611.1180
fval=
3.6576
exitflag=
1
§3二次规划模型
数学模型:
其中H为二次型矩阵,A、Aeq分别为不等式约束与等式约束系数矩阵,f,b,beq,lb,ub,x为向量。
求解二次规划问题函数为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,lb,ub,x0,options)
[x,fval]=quadprog(…)
[x,fval,exitflag]=quadprog(…)
[x,fval,exitflag,output]=quadprog(…)
[x,fval,exitflag,output,lambda]=quadprog(…)
说明:
输入参数中,x0为初始点;若无等式约束或无不等式约束,就将相应的矩阵和向量设置为空;options为指定优化参数。
输出参数中,x是返回最优解;fval是返回解所对应的目标函数值;exitflag是描述搜索是否收敛;output是返回包含优化信息的结构。
Lambda是返回解x入包含拉格朗日乘子的参数。
例1:
求解:
二次规划问题
minf(x)=x1-3x2+3x12+4x22-2x1x2
s.t2x1+x2≤2
-x1+4x2≤3
程序:
f=[1;-3]
H=[6-2;-28]
A=[21;-14]
b=[2;3]
[X,fval,exitflag]=quadprog(H,f,A,b)
结果:
X=
-0.0455
0.3636
fval=
-0.5682
exitflag=
1
例2:
求解:
二次规划问题
min+x12+2x22-2x1x2-4x1-12x2
s.tx1+x2≤2
-x1+2x2≤2
2x1+x2≤3
0≤x1,0≤x2
程序:
H=[2-2;-24];
f=[-4;-12];
A=[11;-12;21];
b=[2;2;3];
lb=zeros(2,1);
[x,fval,exitflag]=quadprog(H,f,A,b,[],[],lb)
结果:
x=
0.6667
1.3333
fval=
-16.4444
exitflag=
1
§4多目标规划模型
多目标规划定义为在一组约束下,多个不同的目标函数进行优化设计。
数学模型:
s.tgj(x)≤0j=1,2,…,k
其中x=(x1,x2,…,xn)为一个n维向量;fi(x)为目标函数,i=1,2,…,m;gj(x)为系统约束,j=1,2,…,k。
当目标函数处于冲突状态时,不存在最优解使所有目标函数同时达到最优。
于是我们寻求有效解(又称非劣解或非支配解或帕累托解)
定义:
若
(
∈Ω)的邻域内不存在Δx,使得(
+Δx∈Ω),且
则称
为有效解。
多目标规划问题的几种常用解法:
(1)
(1) 主要目标法
其基本思想是:
在多目标问题中,根据问题的实际情况,确定一个目标为主要目标,而把其余目标作为次要目标,并且根据经验,选取一定的界限值。
这样就可以把次要目标作为约束来处理,于是就将原来的多目标问题转化为一个在新的约束下的单目标最优化问题。
(2)
(2) 线性加权和法
其基本思想是:
按照多目标fi(x)(i=1,2,…,m)的重要程度,分别乘以一组权系数λj(j=1,2,…,m)然后相加作为目标函数而构成单目标规划问题。
即
,其中
例1:
某钢铁厂准备用5000万用于A、B两个项目的技术改造投资。
设x1、x2分别表示分配给项目A、B的投资。
据专家预估计,投资项目A、B的年收益分别为70%和66%。
同时,投资后总的风险损失将随着总投资和单项投资的增加而增加,已知总的风险损失为0.02x12+0.01x22+0.04(x1+x2)2,问应如何分配资金才能使期望的收益最大,同时使风险损失为最小。
建立数学模型
maxf1(x)=70x1+66x2
minf2(x)=0.02x12+0.01x22+0.04(x1+x2)2
s.tx1+x2≤5000
0≤x1,0≤x2
线性加权构造目标函数:
maxf=0.5f1(x)–0.5f2(x)
化最小值问题:
min(-f)=-0.5f1(x)+0.5f2(x)
首先编辑目标函数M文件ff11.m
functionf=ff11(x)
f=-0.5*(70*x
(1)+66*x
(2))+0.5*(0.02*x
(1)^2+0.01*x
(2)^2+0.04*(x
(1)+x
(2))^2);
调用单目标规划求最小值问题的函数
x0=[1000,1000]
A=[11];
b=5000;
lb=zeros(2,1);
[x,fval,exitflag]=fmincon(@ff11,x0,A,b,[],[],lb,[])
f1=70*x
(1)+66*x
(2)
f2=0.02*x
(1)^2+0.01*x
(2)^2+0.04*(x
(1)+x
(2))^2
结果:
x=
307.1428414.2857
fval=
-1.2211e+004
exitflag=
1
f1=4.8843e+004
f2=2.4421e+004
(3)(3) 极大极小法
其基本思想是:
对于极小化的多目标规划,让其中最大的目标函数值尽可能地小为此,对每个x∈R,我们先求诸目标函数值fi(x)的最大值,然后再求这些最大值中的最小值。
即构造单目标规划:
(4)(4) 目标达到法
对于多目标规划:
s.tgj(x)≤0j=1,2,…,n
先设计与目标函数相应的一组目标值理想化向量
,
再设γ为一松弛因子标量。
设
为权值系数向量。
于是多目标规划问题化为:
在Matlab的优化工具箱中,fgoalattain函数用于解决此类问题。
其数学模型形式为:
minγ
F(x)-weight·γ≤goal
c(x)≤0
ceq(x)=0
Ax≤b
Aeqx=beq
lb≤x≤ub
其中,x,weight,goal,b,beq,lb和ub为向量,A和Aeq为矩阵,c(x),ceq(x)和F(x)为函数,
调用格式:
x=fgoalattain(F,x0,goal,weight)
x=fgoalattain(F,x0,goal,weight,A,b)
x=fgoalattain(F,x0,goal,weight,A,b,Aeq,beq)
x=fgoalattain(F,x0,goal,weight,A,b,Aeq,beq,lb,ub)
x=fgoalattain(F,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)
x=fgoalattain(F,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options)
x=fgoalattain(F,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2)
[x,fval]=fgoalattain(…)
[x,fval,attainfactor]=fgoalattain(…)
[x,fval