MATLAB线性规划(LP)
MATLAB线性规划(LP)
——线性规划的一般形式:目标函数和约束条件都是设计变量的线性函数
1、求解思路
- 结合题目得到目标函数(一般为min型,max型利用求相反值转化为min型)
- 列出约束条件
- 将线性函数函数以及条件转换成矩阵形式
- 编写程序求解
- 要求最优解符合要求
2、例题练习
练习 L 1 \mathscr{L_1} L1
思路
- c数组A矩阵b数组均已知,有 ( x 1 . . . x 6 ) T (x_1...x_6)^T (x1...x6)T;
- 有不等式约束,无等式约束, x ≥ 0 x \geq 0 x≥0.
线性形式
- 目标函数: m i n Z = ∑ i = 1 6 c i x i minZ=\sum_{i=1}^6 c_ix_i minZ=∑i=16cixi
- s . t . = { ∑ k = 1 6 a i k x k ≤ b k , i = 1 , 2... n x i ≥ 0 , i = 1 , 2... n s.t.=\begin{cases} \sum_{k=1}^6a_{ik}x_k \leq b_k ,i=1,2...n\\ x_i \geq 0,i=1,2...n\end{cases} s.t.={ ∑k=16aikxk≤bk,i=1,2...nxi≥0,i=1,2...n
矩阵形式
-
目标函数: m i n Z = c X minZ=cX minZ=cX
-
约束条件: s . t . = { A x ≤ b v l b ≤ x ≤ v u b s.t.=\begin{cases} Ax \leq b \\ vlb \leq x\leq vub\end{cases} s.t.={ Ax≤bvlb≤x≤vub
代码
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)
% minZ=cX z为cX的求和的最小值,c为各项x的系数
% AX≤b,A、X、b均为数组,这是不等式约束条件
% AeqX=beq,Aeq、beq均为数组,为等式约束条件,Aeq为各项x的系数
% 若没有不等式存在,则令A=[],b=[].
% 若没有等式约束 则令Aeq=[], beq=[].
% vlb≤X≤vub,是对变量X的约束
% 返回最优解x以及x处目标函数值fval
答案
x =
1.0e+04 *
3.5000
0.5000
3.0000
0
0
0
fval =
-25000
练习 L 2 \mathscr{L_2} L2
思路
- c数组A矩阵b数组均已知,有 x 1 − x 3 x_1-x_3 x1−x3;
- 有不等式约束,有等式约束, x ≥ 0 x \geq 0 x≥0.
线性形式
- 目标函数: m i n Z = ∑ i = 1 6 c i x i minZ=\sum_{i=1}^6 c_ix_i minZ=∑i=16cixi
- s . t . = { ∑ k = 1 3 a i k x k ≤ b k , i = 1 , 2 , 3 x 1 ≥ 30 x 2 ≥ 0 x 3 ≥ 20 s.t.=\begin{cases} \sum_{k=1}^3a_{ik}x_k \leq b_k ,i=1,2,3\\ x_1 \geq 30 \\ x_2 \geq 0 \\ x_3 \geq 20 \end{cases} s.t.=⎩⎪⎪⎪⎨⎪⎪⎪⎧∑k=13aikxk≤bk,i=1,2,3x1≥30x2≥0x3≥20
矩阵形式
-
目标函数: m i n Z = c X minZ=cX minZ=cX
-
约束条件: s . t . = { A x ≤ b v l b ≤ x ≤ v u b A e q x = b e q s.t.=\begin{cases} Ax \leq b \\ vlb \leq x\leq vub \\ Aeq x = beq \end{cases} s.t.=⎩⎪⎨⎪⎧Ax≤bvlb≤x≤vubAeqx=beq
代码
c=[6 3 4];
A=[0 1 0];
b=[50];
vlb=[30 0 20];
vub=[];
Aeq=[1 1 1];
beq=[120];
[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub);
答案
x =
30
50
40
fval =
490
3、实际建模
3.1 提出问题
3.2 基本假设&符号规定
3.2.1 基本假设
- 总投资数额M很大,为了计算方便将M设为1;
- n种资产 s i s_i si相互独立,互不相干;
- 投资越分散,总体风险越小;
- 银行存款无交易费也无风险( p 0 = q 0 = 0 p_0=q_0=0 p0=q0=0);
- 在投资期间的 r i , q i , p i , r 0 r_i,q_i,p_i,r_0 ri,qi,pi,r0不受外界环境因素干扰,均为定值;
- 净收益与总体风险只与 r i , q i , p i , r 0 r_i,q_i,p_i,r_0 ri,qi,pi,r0相关,不受其他无关因素干扰;
- 总体风险用n种投资项目 S i S_i Si中风险最大的一个来度量。
3.2.2 符号规定
- S i S_i Si —— 投资的第 i i i种资产
- M M M —— 投资总额
- r i q i p i r_i\ q_i\ p_i ri qi pi —— 某一资产的平均收益率、风险损失率、交易费率
- r 0 r_0 r0 —— 同期银行的存款利率
- u i u_i ui —— 某一资产 S i S_i Si的交易定额
- x i x_i xi —— 某一资产 S i S_i Si的投资额
- a a a —— 投资风险度
- k k k —— 可接受的最低总收益
- Q Q Q —— 总体收益
- Δ Q \Delta Q ΔQ—— 总体收益的增量
- L L L —— 总体风险损失
- s s s —— 权重因素
3.3 建立并分析模型
3.3.1 基本模型
-
总体风险用n种投资项目 S i S_i Si中风险最大的一个来度量,则 m a x { q i x i ∣ i = 1 , 2... n } max\{q_ix_i|i=1,2...n\} max{ qixi∣i=1,2...n};
-
购买 S i Si Si( i = 1 , 2... n i=1,2...n i=1,2...n)的交易费为:
{ p i x i , x i ≥ u i p i u i , x i < u i \begin{cases} p_ix_i,x_i \geq u_i \\ p_iu_i,x_i<u_i\end{cases} { pixi,xi≥uipiui,xi<ui
且投资数额 M M M相较于 p i u i p_iu_i piui非常大,则可忽略不计,交易费记为 p i x i p_ix_i pixi,则购买某资产 S i S_i Si的净收益为 x i ( r i − p i ) x_i(r_i-p_i) xi(ri−pi); -
投资越分散,总体风险越小,所以总体收益的线性规划是一个多目标的线性规划,
目标函数:
Q = m a x ∑ i = 0 n x i ( r i − p i ) L = m i n { m a x { q i x i } } , ( i = 1 , 2... n ) Q=max{\sum_{i=0}^nx_i(r_i-p_i)} \\ L=min\{max\{q_ix_i\}\},(i=1,2...n) Q=maxi=0∑nxi(ri−pi)L=min{ max{ qixi}},(i=1,2...n)
ps:根据不同的建模需求,目标函数会被一定程度地组合成新的目标函数,或被改造成约束条件。
约束条件:
∑ i = 0 n x i ( 1 + p i ) = M x i ≥ 0 , ( i = 0 , 1... n ) \sum_{i=0}^{n}x_i(1+p_i)=M \\ x_i \geq 0,(i=0,1...n) i=0∑nxi(1+pi)=Mxi≥0,(i=0,1...n)
3.3.2 简化模型
模型1
固定风险水平,优化收益
要求:
- 满足总体风险度 max { q i x i } / M ≤ a , ( i = 1 , 2... n ) \max\{q_ix_i\}/M \leq a,(i=1,2...n) max{ qixi}/M≤a,(i=1,2...n),且作为限制条件;
- 目标函数为Q.
模型2
固定收益,降低风险
要求:
- 满足总体收益 Q = ∑ i = 0 n x i ( r i − p i ) ≥ k Q={\sum_{i=0}^nx_i(r_i-p_i)} \geq k Q=∑i=0nxi(ri−pi)≥k,且作为限制条件;
- 目标函数为L.
模型3
根据赋予的权重因素s进行预期收益与资产风险的取舍
要求:
- 约束条件同基本模型
- 0 < s ≤ 1 0< s \leq 1 0<s≤1
- 此题中 s s s为投资偏好系数,目标函数为 s L − ( 1 − s ) Q sL-(1-s)Q sL−(1−s)Q.
3.4 求解模型
以模型1为例:
-
目标函数: m i n f v a l = ( − 0.05 − 0.27 − 0.19 − 0.185 − 0.185 ) ( x 0 x 1 x 2 x 3 x 4 ) T minfval = (-0.05\ -0.27\ -0.19\ -0.185\ -0.185)(x_0\ x_1\ x_2\ x_3\ x_4)^T minfval=(−0.05 −0.27 −0.19 −0.185 −0.185)(x0 x1 x2 x3 x4)T
-
约束条件: { x 0 + 0.01 x 1 + 0.02 x 2 + 0.045 x 3 + 0.065 x 4 = M = 1 0.025 x 1 ≤ a 0.015 x 2 ≤ a 0.055 x 3 ≤ a 0.026 x 4 ≤ a x i ≥ 0 , ( i = 0 , 1...4 ) \begin{cases} x_0+0.01x_1+0.02x_2+0.045x_3+0.065x_4=M=1 \\ 0.025x_1 \leq a \\ 0.015x_2 \leq a \\ 0.055x_3 \leq a \\ 0.026x_4 \leq a \\ x_i \geq 0,(i=0,1...4) \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧x0+0.01x1+0.02x2+0.045x3+0.065x4=M=10.025x1≤a0.015x2≤a0.055x3≤a0.026x4≤axi≥0,(i=0,1...4)
-
a a a不定,则令 a = 0 a=0 a=0,步长 Δ a = 0.0001 \Delta a=0.0001 Δa=0.0001进行循环搜索,得到结果,编写函数如下:
a=0; while (1.1-a)>1 c=[-0.05 -0.27 -0.19 -0.185 -0.185]; A=[0 0.025 0 0 0;0 0 0.015 0 0;0 0 0 0.055 0;0 0 0 0 0.026]; b=[a;a;a;a]; Aeq=[1 1.01 1.02 1.045 1.065]; beq=[1]; vlb=[0;0;0;0;0]; vub=[]; [x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub); Q=-fval; % 负值才为净利益 plot(a,Q,'-r.'); % 制图 axis([0 0.1 0 0.5]);% 坐标轴 hold on; % 不清理前面循环过程的作图 a=a+0.001; % 循环 end xlabel('a'); ylabel('Q'); title('模型1求解');
-
结果
3.5 总结与分析
-
风险越大,利益越大;
-
分析图中散点图可知,投资越分散,风险越小,投资越密集,风险越大;
-
(a,Q)=(0.006,0.2019)为图中曲线的拐点,拐点往左,风险小利益较少,但在持续增加,拐点往右,风险大,且利益增长缓慢,如果对利益没有特殊追求追求的话,应该选择该点作为投资参照;
-
则模型1的最优解为:
x =
0
0.2400
0.4000
0.1091
0.2212
Q =
0.2019
4、实验作业
4.1 提出问题
4.2 建立模型
4.2.1 思路
- Q Q Q单位万元, x i x_i xi单位百箱;
- 使用min模型求解;
- 有不等式约束,无等式约束,x均为正数。
4.2.2 线性形式
- 目标函数: m a x Q = 10 x 1 + 9 x 2 maxQ=10x_1+9x_2 maxQ=10x1+9x2
- s . t . = { 6 x 1 + 5 x 2 ≤ 60 10 x 1 + 20 x 2 ≤ 150 0 ≤ x 1 ≤ 8 0 ≤ x 2 s.t.=\begin{cases} 6x_1+5x_2 \leq 60 \\ 10x_1+20x_2 \leq 150 \\ 0 \leq x_1 \leq 8 \\ 0 \leq x_2 \end{cases} s.t.=⎩⎪⎪⎪⎨⎪⎪⎪⎧6x1+5x2≤6010x1+20x2≤1500≤x1≤80≤x2
4.2.3 矩阵形式
-
目标函数: m a x Q = ∑ i = 1 2 c X maxQ=\sum_{i=1}^{2}cX maxQ=∑i=12cX
-
约束条件: s . t . = { A x ≤ b v l b ≤ x ≤ v u b s.t.=\begin{cases} Ax \leq b \\ vlb \leq x\leq vub\end{cases} s.t.={ Ax≤bvlb≤x≤vub
4.3 求解模型
4.3.1 代码
% c=[-10 -9];
% A=[6 5;10 20;1 0];
% b=[60;150;8];
% Aeq=[];
% beq=[];
% vlb=[0;0];
% vub=[];
% [x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub);
% Q=-fval;
a=60
while(62-a>0)
c=[-10 -9];
A=[6 5;10 20;1 0];
b=[a;150;8];
Aeq=[];
beq=[];
vlb=[0;0];
vub=[];
[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub);
Q=-fval;
plot(a,Q,'-r.');% 制图
axis on; % 坐标轴
hold on; % 不清理前面循环过程的作图
a=a+0.01;
end
xlabel('a');
ylabel('Q');
% c=[-11 -9];
% A=[6 5;10 20;1 0];
% b=[60;150;8];
% Aeq=[];
% beq=[];
% vlb=[0;0];
% vub=[];
% [x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub);
% Q=-fval;
4.3.2 答案
问题0:考虑到工人人数为整数,所以应该左右寻找满足 10 x 1 10x_1 10x1以及 5 x 2 5x_2 5x2为整数的 x 1 、 x 2 x_1、x_2 x1、x2.
x =
6.4286
4.2857
Q =
102.8571
%优化后
x =
6.5000
4.2000
Q =
102.8000
问题1:由图知,应该作这项投资。
问题2:修改系数可得 x 、 Q x、Q x、Q,且正好满足 10 x 1 10x_1 10x1以及 5 x 2 5x_2 5x2为整数,则应该改变生产计划。
x =
8.0000
2.4000
Q =
109.6000