MATLAB线性规划

demo1

clc, clear
c = [4;3]; b = [10;8;7];
a = [2,1;1,1;0,1]; lb = zeros(2,1);
[x,fval] = linprog(-c,a,b,[],[],lb) %没有等号约束
y = -fval  %目标函数为最大化
1
2
3
4
5

clc, clear
prob = optimproblem('ObjectiveSense', 'max');
%这里ObjectiveSense是目标类型为求最大值的意思
%有时候直接是prob=optimproblem就是默认求最小值

c = [4;3]; b = [10;8;7];
a = [2,1;1,1;0,1]; lb = zeros(2,1);
x = optimvar('x',2,'LowerBound',0);

%x=optimvar('x',1,2,'TYPE','integer','LowerBound',0,'UpperBound',inf);
%optimvar函数是一种类似于赋值的函数
%第一个‘x’里面是变量名，后面说的是该变量所包含的行数和列数
%‘TYPE’，后面定义的是该函数所属类型，比如说integer整数型，double双精度型号等
%‘LowerBound'与'UpperBound'表示下界与上界所跟的0，inf分别是范围

prob.Objective = c'*x;
prob.Constraints.con = a*x<=b;
[sol, fval, flage, out] = solve(prob)
sol.x  %显示决策变量的值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

fval是得到规划后的y值，sol是储存最优x解的一个元胞

sol，fval，flag，out】=solve（prob）
xx=sol.x，yy=sol.y
【i，j】=find（yy）；i j=【i；j】

参考博文

demo2

clc, clear, f=[-2; -3; 5];
a=[-2,5,-1;1,3,1]; b=[-10;12];
aeq=[1,1,1]; beq=7;
[x,y]=linprog(f,a,b,aeq,beq,zeros(3,1));
x, y=-y  %目标函数最大化

1
2
3
4
5
6

clc, clear
prob = optimproblem('ObjectiveSense','max');
x = optimvar('x',3,'LowerBound',0);
prob.Objective = 2*x(1) + 3*x(2)-5*x(3);
prob.Constraints.con1 = x(1)+x(2)+x(3) == 7;
prob.Constraints.con2 = 2*x(1)-5*x(2)+x(3) >=10;
prob.Constraints.con3 = x(1)+3*x(2)+x(3) <=12;
[sol,fval,flag,out]= solve(prob), sol.x

1
2
3
4
5
6
7
8
9

demo3

clc, clear
prob = optimproblem
x = optimvar('x',4,4,'LowerBound',0);
prob.Objective = 2800*sum(x(:,1))+4500*sum(x(1:3,2))+...
    6000*sum(x(1:2,3))+7300*x(1,4);
prob.Constraints.con = [sum(x(1,:))>=15,
    sum(x(1,2:4))+sum(x(2,1:3))>=10,
    x(1,3)+x(1,4)+x(2,2)+x(2,3)+x(3,1)+x(3,2)>=20,
    x(1,4)+x(2,3)+x(3,2)+x(4,1)>=12];
[sol,fval,flag,out]= solve(prob), sol.x

1
2
3
4
5
6
7
8
9
10
11

投资组合模型(ex1_9_1)

clc, clear, close all
prob = optimproblem('ObjectiveSense','max');
x = optimvar('x',5,1,'LowerBound',0);
c=[0.05,0.27,0.19,0.185,0.185]; %净收益率
Aeq=[1,1.01,1.02,1.045,1.065];  %等号约束矩阵
prob.Objective = c*x; M = 10000; 
prob.Constraints.con1 = Aeq*x==M; %等号约束条件
q=[0.025,0.015,0.055,0.026]'; %风险损失率
a = 0; aa = []; QQ = []; XX = []; hold on
while a<0.05
    prob.Constraints.con2 = q.*x(2:end)<=a*M;
    [sol,Q,flag,out]= solve(prob);
    aa = [aa; a]; QQ = [QQ,Q]; 
    XX = [XX; sol.x']; a=a+0.001;
end
plot(aa, QQ, '*k')
xlabel('$a$','Interpreter','Latex'),
ylabel('$Q$','Interpreter','Latex','Rotation',0)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

clc, clear, close all, format long g
M =10000; prob = optimproblem; 
x = optimvar('x',6,1,'LowerBound',0);
r=[0.05,0.28,0.21,0.23,0.25]; %收益率
p=[0, 0.01, 0.02, 0.045, 0.065]; %交易费率
q=[0, 0.025, 0.015, 0.055, 0.026]'; %风险损失率
%w = 0:0.1:1 
w = [0.766, 0.767, 0.810, 0.811, 0.824, 0.825, 0.962, 0.963, 1.0]
V = []; %风险初始化
Q = []; %收益初值化
X = []; %最优解的初始化
prob.Constraints.con1 = (1+p)*x(1:end-1)==M;
prob.Constraints.con2 = q(2:end).*x(2:end-1)<=x(end);
for i = 1:length(w)
    prob.Objective = w(i)*x(end)-(1-w(i))*(r-p)*x(1:end-1);
    [sol,fval,flag,out]=solve(prob);
    xx = sol.x; V=[V,max(q.*xx(1:end-1))];
    Q=[Q,(r-p)*xx(1:end-1)]; X=[X; xx'];
    plot(V, Q, '*-'); grid on
    xlabel('风险（元）'); ylabel('收益（元）')
end
V, Q, format


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

demo

clc, clear, prob = optimproblem;
x = optimvar('x',2,'LowerBound',0);
prob.Objective = sum(x.^2)-4*x(1)+4;
con = [-x(1)+x(2)-2 <= 0
x(1)^2-x(2)+1 <= 0];  %不等式约束
prob.Constraints.con = con;
x0.x = rand(2,1)  %非线性规划必须赋初值
[sol,fval,flag,out] = solve(prob,x0), sol.x

1
2
3
4
5
6
7
8
9

x0.x = rand(2,1)  %非线性规划必须赋初值
1

简单非线性规划模型之 彩电销售

clc, clear, format long g
syms x1 x2  %定义符号变量
f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2)
f=simplify(f) %化简目标函数
f1=diff(f,x1), f2=diff(f,x2)  %求目标函数关于x1,x2的偏导数
[x10,x20]=solve(f1,f2)  %解代数方程求驻点
x10=round(double(x10))   %取整
x20=round(double(x20))   %取整
f0=subs(f,{x1,x2},{x10,x20})  %求目标函数的取值
f0=double(f0)
subplot(121), fmesh(f,[0,10000,0,10000]),title('') %画三维图形
xlabel('$x_1$','Interpreter','Latex')
ylabel('$x_2$','Interpreter','Latex')
subplot(122), c=fcontour(f,[0,10000,0,10000]) %绘制等高线
contour(c.XData, c.YData, c.ZData, 'ShowText','on') %等高线标注
xlabel('$x_1$','Interpreter','Latex')
ylabel('$x_2$','Interpreter','Latex')
p1=339-0.01*x10-0.003*x20  %计算19英寸的平均售价
p2=399-0.004*x10-0.01*x20  %计算21英寸的平均售价
c=400000+195*x10+225*x20  %计算总支出
rate=f0/c  %计算利润率
% 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23