数学建模算法与应用（司守奎）Matlab版第一章线性规划代码解析

最新推荐文章于 2025-05-04 10:19:44 发布

你干码，哎哟

最新推荐文章于 2025-05-04 10:19:44 发布

阅读量687

点赞数 18

文章标签：数学建模算法 matlab

本文链接：https://blog.youkuaiyun.com/cyc520521/article/details/145145948

版权

代码仅给出了部分注释，详细的问题解决方法请查阅《数学建模算法与应用》以及《数学建模算法与应用习题解答》，有不懂的地方欢迎在评论区提问，看到必答！

1.1

1.

clc,clear
f=[-4,-3];%求最大值将f变号，以及结果的值变号
A=[2,1;1,1;0,1];b=[10,8,7];
lb=zeros(2,1);%两行一列的零矩阵
[x,y]=linprog(f,A,b,[],[],lb);%无等号约束用空矩阵
x,W=-y %输出结果

2.

clc,clear
prob=optimproblem('ObjectiveSense','max');%创建目标函数最大化的优化问题
f=[4,3];A=[2,1;1,1;0,1];b=[10;8;7];%换成了列向量
x=optimvar('x',2,'LowerBound',0);%决策变量
prob.Objective = f*x;%目标函数，矩阵乘法
prob.Constraints.con=A*x<=b;%约束条件
[sol,fval,flag,out]=solve(prob);%fval返回最优值
sol.x,fval, %显示结果

1.3

f=[-2,-3,5];
A=[-2,5,-1;1,3,1];b=[-10,12];%每一个不等式的系数是一行
Aeq=ones(1,3);Beq=7;
lb=zeros(3,1);
[x,y]=linprog(f,A,b,Aeq,Beq,lb);
x,y=-y

1.4

clc,clear
prob=optimproblem;%('ObjectiveSense','min');%默认最小化
x=optimvar('x',4,4,'LowerBound',0);%4*4个x
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,1:4))>=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,fval

1.5

txt文件

6,2,6,7,4,2,5,9,60
4,9,5,3,8,5,8,2,55
5,2,1,9,7,4,3,3,51
7,6,7,3,9,2,7,1,43
2,3,9,5,7,2,6,5,41
5,5,2,2,8,1,4,3,52
35,37,22,32,41,32,43,38,0

clc,clear,a=load('wenti.txt');%将表中的信息放在文本中
c=a(1:end-1,1:end-1);%（行取1到倒数第2，列取1到倒数第2）
e=a(1:end-1,end);d=a(end,1:end-1);%e为产量，d为需求量
prob=optimproblem;%默认最小化
x=optimvar('x',6,8,'LowerBound',0);%决策变量x6*8个，下限为0
prob.Objective=sum(sum(c.*x));%目标函数：运价*运送量再求和
%里面的sum求每列的和返回一行值，外面的sum求这一行的和
prob.Constraints.con1=sum(x,1)==d;%第一维按列求和等于需求
prob.Constraints.con2=sum(x,2)<=e;%第二维按行求和大于等于产量
[sol,fval,flag,out]=solve(prob);
fval,xx=sol.x,%显示结果
writematrix(xx,'answer.xlsx')%将xx的数值放到文本中

1.7

clc,clear
c=[1:4];b=[-2;-1;-1/2];
a=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3];
prob=optimproblem;
u=optimvar('u',4,'LowerBound',0);
v=optimvar('v',4,'LowerBound',0);
prob.Objective=c*(u+v);%矩阵乘法
prob.Constraints.con=a*(u-v)<=b;
[sol,fval,flag,out]=solve(prob);
fval,x=sol.u+sol.v

1.9

%模型一
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];%1+费用率
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%使图像叠加而不是消除,hoid off取消
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];%;将a添加到此列下方50*1,,将Q添加到此行右方1*50
    XX=[XX;sol.x'];a=a+0.001;
end
plot(aa,QQ,'*k')%二维绘图,plot(横坐标,纵坐标，星号标记黑色black)
xlabel('$a$','Interpreter','Latex')%给坐标轴设置名称，Interpreter用于指定Latex解释器
ylabel('$Q$','Interpreter','Latex','Rotation',0)%Rotation，0，防止y轴名称翻转

%模型二
clc,clear,close all,format long g%长浮点数表示，结尾format恢复默认
M=10000;prob=optimproblem;%求最小值
x=optimvar('x',6,1,'LowerBound',0);%x(6)为风险损失的最大值
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);%x(6)为最大值
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%显示坐标网格线
end
xlabel('风险/元');ylabel('收益/元')
V,Q,format

%模型三
clc,clear,close all
prob=optimproblem('ObjectiveSense','min');
c=[0.05,0.27,0.19,0.185,0.185];
q=[0.025,0.015,0.055,0.026];
Aeq=[1,1.01,1.02,1.045,1.065];
x=optimvar('x',6,1,'LowerBound',0);
prob.Objective=x(6);M=10000;
prob.Constraints.con1=Aeq*x(1:end-1)==M;
prob.Constraints.con2=q'.*x(2:end-1)<=x(6);
k=0.05;kk=[];QQ=[];XX=[];hold on
while k<0.26
    prob.Constraints.con3=c*x(1:end-1)>=k*M;
    [sol,fval,flag,out]=solve(prob);
    kk=[kk;k];QQ=[QQ;fval];XX=[XX;sol.x];
    k=k+0.005;
end
plot(kk,QQ,"*-k")
xlabel('$k$','Interpreter','latex');
ylabel('风险/元');

T1

clc,clear
%一
f=[3,-1,-1];
A=[1,-2,1;4,-1,-2];b=[11;-3];
Aeq=[-2,1];
prob=optimproblem('ObjectiveSense','max');
x=optimvar('x',3,'LowerBound',0);
prob.Objective=f*x;
prob.Constraints.con1=A*x<=b;
prob.Constraints.con2=Aeq*[x(1);x(3)]==1;
[sol,fval,flag,out]=solve(prob);
fval,sol.x

clc,clear
%二
f=[-3,1,1];
A=[1,-2,1;4,-1,-2];b=[11;-3];
Aeq=[-2,0,1];beq=1;
lb=zeros(3,1);
[x,y]=linprog(f,A,b,Aeq,beq,lb);
fval=-y,x

T2

clc,clear
c=[1:4];
Aeq=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3];beq=[0,1,-1/2];
prob=optimproblem;
u=optimvar('u',4,'LowerBound',0);
v=optimvar('v',4,'LowerBound',0);
prob.Objective=c*(u+v);
prob.Constraints.con=Aeq*(u-v)==beq';
[sol,fval,flag,out]=solve(prob);
fval,x=sol.u-sol.v

T3

clc,clear,format long g%整数规划
prob=optimproblem("ObjectiveSense","max");
x = optimvar('x',9,'Type','integer','LowerBound',0);%所有解是整数，x(10)=x(9)
prob.Objective=(x(1)+x(2))+x(8)*1.65+x(9)*2.3-...%换行用...
(5*x(1)+10*x(6))*0.05-(7*x(2)+9*x(7)+12*x(9))*0.0321-...
(6*x(3)+8*x(8))*250/4000-(4*x(4)+11*x(9))*783/7000-x(5)*0.35;
prob.Constraints.con1=[x(1)+x(2)==x(3)+x(4)+x(5)
    x(6)+x(7)==x(8)];%每种刚好够
prob.Constraints.con2=[5*x(1)+10*x(6)<=6000
    7*x(2)+9*x(7)+12*x(9)<=10000
    6*x(3)+8*x(8)<=4000
    4*x(4)+11*x(9)<=7000
    7*x(5)<=4000];%不能超过机器有效时长
[sol,fval,flag]=solve(prob);
x=fval,sol.x,format

T4

clc,clear
%货物i送到舱j中i=1,2,3,4;j=1,2,3
c=[3100,3800,3500,2850];%利润
mlimt1=[18,15,23,12];v=[480,650,580,390];
mlimt2=[10,16,8];vlimt=[6800,8700,5300];
prob=optimproblem("ObjectiveSense","max");
x=optimvar('x',4,3,'LowerBound',0);
A=sum(x,1);%货舱所放货物的吨数
B=sum(x,2);%运送每个货物的吨数
prob.Objective=B'*c';%总利润
prob.Constraints.con1=B'<=mlimt1;%货物吨数限制
prob.Constraints.con2=A<=mlimt2;%货舱吨数限制
prob.Constraints.con3=v*x<=vlimt;%货舱体积限制
prob.Constraints.con4=[%质量比例限制
    A(1)/5==A(2)/8
    A(2)/2==A(3)];
[sol,fval,flag,out]=solve(prob);
fval,x=sol.x

T5

clc,clear,format long g
prob=optimproblem('ObjectiveSense','max');
x=optimvar('x',4,5,'LowerBound',0);
prob.Objective=1.15*x(1,4)+1.25*x(2,3)+1.4*x(3,2)+1.06*x(4,5);
con1=[x(1,1)+x(4,1)==100000 %最初十万第一年投完
    x(1,2)+x(3,2)+x(4,2)==1.06*x(4,1) %后一年投资等于前一年剩余
    x(1,3)+x(2,3)+x(4,3)==1.15*x(1,1)+1.06*x(4,2)
    x(1,4)+x(4,4)==1.15*x(1,2)+1.06*x(4,3)
    x(4,5)==1.15*x(1,3)+1.06*x(4,4)];
con2=[x(2,3)<=40000
    x(3,2)<=30000];%条件约束
prob.Constraints.con1=con1;
prob.Constraints.con2=con2;
[sol,fval,flag,out]=solve(prob);
fval,x=sol.x,format

T6

clc,clear
p=[24,15];c=[20,12,8];
mlimt1=[600,750,625];mlimt2=[600,800];
prob=optimproblem("ObjectiveSense","max");
x=optimvar('x',2,3,'LowerBound',0);
A=sum(x,1);B=sum(x,2);
prob.Objective=sum(B.*p')-sum(A.*c);
prob.Constraints.con1=[%别用除法
    B(1)*0.5<=x(1,1)
    B(1)*0.25<=x(1,2)
    x(1,3)<=0.1*B(1)
    x(2,1)<=0.4*B(2)
    x(2,2)<=0.4*B(2)
    B(2)*0.15<=x(2,3)];
prob.Constraints.con2=A(1:3)<=mlimt1;%原料限制
prob.Constraints.con3=B(1:2)>=mlimt2';%需求量限制
[sol,fval]=solve(prob);
fval,x=sol.x

T7

clc,clear
rng(0)%进行一致性比较，生成相同的随机数
a=randi([0,10],100,150);%a为100*150，0到10之间的随机整数
prob=optimproblem('ObjectiveSense','max');
v=optimvar('v');x=optimvar('x',100,'LowerBound',0);
prob.Objective=v;
prob.Constraints.con1=a'*x>=v;
prob.Constraints.con2=sum(x)==100;
[sol,fval]=solve(prob);v=sol.v