这篇文章主要介绍了python数学建模基础教程，具有一定借鉴价值，需要的朋友可以参考下。希望大家阅读完这篇文章后大有收获，下面让小编带着大家一起了解一下。
🚧 线性规划
声明：
-
采用如下求解方法：
- python cxpy
- matlab
- lingo
-
5.1 求解下列线性规划的解：
m a x z = 8 x 1 − 2 x 2 + 3 x 3 − 4 x 4 − 2 x 5 { x 1 + x 2 + x 3 + x 4 + x 5 ≤ 400 , x 1 + 2 x 2 + 2 x 3 + x 4 + 6 x 5 ≤ 800 , 2 x 1 + x 2 + 6 x 3 ≤ 200 , x 3 + x 4 + 5 x 5 ≤ 200 , 0 ≤ x i ≤ 99 , i = 1 , 2 , 3 , 4 x 5 ≥ 10 max~z=8x_1-2x_2+3x_3-4x_4-2x_5\\ \left\{ \begin{aligned} x_1+x_2+x_3+x_4+x_5\le&400,\\ x_1+2x_2+2x_3+x_4+6x_5&\le800,\\ 2x_1+x_2+6x_3\le200,\\ x_3+x_4+5x_5\le 200,\\ 0\le x_i \le99,~i=1,2,3,4\\ x_5\ge10 \end{aligned} \right. max z=8x1−2x2+3x3−4x4−2x5⎩ ⎨ ⎧x1+x2+x3+x4+x5≤x1+2x2+2x3+x4+6x52x1+x2+6x3≤200,x3+x4+5x5≤200,0≤xi≤99, i=1,2,3,4x5≥10400,≤800,- matlab
c=-[8,-2,3,-1,-2];%价值向量 A=[ones(1,5);1,2,2,1,6;2,1,6,0,0;0,0,1,1,5]; b=[400,800,200,200]'; [x,fval]=linprog(c,A,b,[],[],[zeros(1,4),-10]',[ones(1,4)*99,inf]'); disp('最优解为：');disp(x) disp('最大值为：');disp(-fval) %结果 Optimal solution found. 最优解为： 99.0000 0 0.3333 0 -10.0000 最大值为： 813- python
#导入凸优化库 import cvxpy as cv import numpy as np #构造系数矩阵 A=np.array([[1]*5,[1,2,2,1,6],[2,1,6,0,0],[0,0,1,1,5]]) b=np.array([400,800,200,200]) c=np.array([8,-2,3,-1,-2]) x=cv.Variable(5)#5个决策变量 x等同与列向量 obj=cv.Maximize(c@x)#构造目标函数 con=[A@x<=b,x[0:-1]<=99*np.ones(4),x[4]>=-10,x[0:4]>=0] prob=cv.Problem(obj,con) prob.solve(solver='GLPK_MI') print('可行解=\n') for i in x.value: print(i,'\n') print('最大值:',prob.value) #结果： 可行解= 99.0 -0.0 0.3333333333333333 -0.0 -10.0 最大值: 813.0- lingo
结果： 813
sets: factory/1..5/:x,c; plant/1..4/:k; endsets data: c=8,-2,3,-1,-2; enddata max=@sum(factory(i):x(i)*c(i)); @sum(factory(i):x(i))<400; x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)<800; 2*x(1)+x(2)+6*x(3)<200; x(3)+x(4)+5*x(5)<200; @for(plant(i):@bnd(0,x(i),99)); @free(x(5)); x(5)>-10; X( 1) 99.00000 X( 2) 0.000000 X( 3) 0.3333333 X( 4) 0.000000 X( 5) -10.00000 -
5.2 求5.4节模型三的解
模型三：
m i n s { m a x 1 ≤ i ≤ n { q i x i } } − ( 1 − s ) ∑ i = 0 n ( r i − p i ) x i { ∑ i = 0 n ( 1 + p i ) x i = M , x i ≥ 0 , i = 0 , 1 , 2 , … min~s\big\{\mathop{max}\limits_{1\le i\le n}\{q_ix_i\}\big\}-(1-s)\sum_{i=0}^{n}(r_i-p_i)x_i\\ \left\{ \begin{aligned} &\sum_{i=0}^{n}(1+p_i)x_i=M,\\ &x_i\ge0,i=0,1,2,\dots \end{aligned} \right. min s{1≤i≤nmax{qixi}}−(1−s)i=0∑n(ri−pi)xi⎩ ⎨ ⎧i=0∑n(1+pi)xi=M,xi≥0,i=0,1,2,…
令 { m a x 1 ≤ i ≤ n { q i x i } } = x 5 \big\{\mathop{max}\limits_{1\le i\le n}\{q_ix_i\}\big\}=x_5 {1≤i≤nmax{qixi}}=x5 ：
m i n f = s x 5 − ( 1 − s ) [ 0.05 , 0.27 , 0.19 , 0.185 , 0.185 ] [ x 0 , x 1 , x 2 , x 3 , x 4 ] T { 0.025 x 1 − x 5 ≤ 0 0.015 x 2 − x 5 ≤ 0 0.055 x 3 − x 5 ≤ 0 0.026 x 4 − x 5 ≤ 0 x 0 + 1.01 x 2 + 1.02 x 2 + 1.045 x 3 + 1.065 x 4 = 1 min \quad f= sx_5-(1-s)[0.05,0.27,0.19,0.185,0.185][x_0,x_1,x_2,x_3,x_4]^T\\ \left\{\begin{aligned} &0.025x_1-x_5\le0\\ &0.015x_2-x_5\le0\\ &0.055x_3-x_5\le0\\ &0.026x_4-x_5\le0\\ &x_0+1.01x_2+1.02x_2+1.045x_3+1.065x_4=1 \end{aligned} \right. minf=sx5−(1−s)[0.05,0.27,0.19,0.185,0.185][x0,x1,x2,x3,x4]T⎩ ⎨ ⎧0.025x1−x5≤00.015x2−x5≤00.055x3−x5≤00.026x4−x5≤0x0+1.01x2+1.02x2+1.045x3+1.065x4=1- python
#利用cxpy import cvxpy as cp import numpy as np import matplotlib.pyplot as plt s = 0.01 point = [] S=[] v = np.array([0.05, 0.27, 0.19, 0.185, 0.185]).reshape(1, 5) X = cp.Variable((6, 1)) # 决策变量 列向量 q = np.array([0.025, 0.015, 0.055, 0.026]) Aeq = np.array([1, 1.01, 1.02, 1.045, 1.065, 0]).reshape(1, 6) # 约束条件 con = [Aeq@X == 1, cp.multiply(q, X[1:5, 0]) - X[5, 0] <= 0,X>=np.zeros((6,1))] while s < 1: f = cp.Minimize(s * X[5, 0] - (1 - s) * (v @ X[0:5,0])) # 构造目标函数 prob=cp.Problem(f,con) prob.solve(solver='GLPK_MI') point.append(-prob.value) S.append(s) s+=0.02 plt.plot(S,point,'ro') plt.show()
- matlab
s=0.01; V=[0.025,0.015,0.005,0.026]; %构造系数矩阵 A=diag(V,1);A(end,:)=[]; A(:,6)=-ones(4,1); Aeq=[1,1.01,1.02,1.045,1.065,0]; beq=1; while s<1 %构造价值向量 c=(s-1).*[0.05,0.27,0.19,0.185,0.185,0]; c(6)=s; %求解最优解 [x,fval]=linprog(c,A,zeros(4,1),Aeq,beq,zeros(6,1),inf.*ones(6,1)); s=s+0.025; hold on plot(s,-fval,'ko') hold off end
- lingo 在此处无法实现可视化操作，所以建议使用上述的两种语言
-
5.3 某股民决定对6家公司的股票进行投资，根据对着6家公司的了解，估计这6家公司股票的明年预期收益和这6种股票的收益的协方差矩阵的数据见表5.6.要获得至少25%的预期收益，最小风险是多少
- 目标函数：风险 R a n k ( p ) = w ′ V w , w 为投资组合的权重， V 为协方差矩阵 Rank(p)=w^{'}Vw,w为投资组合的权重，V为协方差矩阵 Rank(p)=w′Vw,w为投资组合的权重，V为协方差矩阵
- 决策变量：权重向量的元素（6个）
- 约束条件：权重之和=1；收益 ≤ 25 % \le 25\% ≤25%
- 二次规划问题
m i n R a n k = [ x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ] T ⋅ V ⋅ [ x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ] { ∑ i = 1 6 x i = 1 , i = 1 , 2 , … 6 0.02 x 1 + 0.42 x 2 + x 3 + 0.5 x 4 + 0.46 x 5 + 0.3 x 6 ≥ 0.25 x i ≥ 0 , i = 1 , 2 , … 6 min~Rank=[x_1,x_2,x_3,x_4,x_5,x_6]^T\cdot V\cdot[x_1,x_2,x_3,x_4,x_5,x_6]\\ \left\{\begin{aligned} & \sum_{i=1}^{6}x_i=1,~~i=1,2,\dots6\\ &0.02x_1+0.42x_2+x_3+0.5x_4+0.46x_5+0.3x_6\ge 0.25\\ &x_i\ge0,~~i=1,2,\dots6 \end{aligned} \right. min Rank=[x1,x2,x3,x4,x5,x6]T⋅V⋅[x1,x2,x3,x4,x5,x6]⎩ ⎨ ⎧i=1∑6xi=1, i=1,2,…60.02x1+0.42x2+x3+0.5x4+0.46x5+0.3x6≥0.25xi≥0, i=1,2,…6
- python
#导入凸优化库 import cvxpy as cv import numpy as np from scipy.optimize import linprog import pandas as pd d=pd.read_excel('协方差矩阵.xlsx') data=np.vstack((np.array(d.columns),d.values)) #六个决策变量 x=cv.Variable((6,1)) #构造目标函数 rank=cv.quad_form(x,data)# """ Alias for :math:`x^T P x`构造矩阵（二次型） #构造约束条件 A=np.array([0.02,0.42,1,0.5,0.46,0.3]).reshape(1,6) con=[cv.sum(x)==1,A@x>=0.25] prob=cv.Problem(cv.Minimize(rank),con) prob.solve() print('可行解为：\n') for i in x.value: print(i,' ') print('最小值:',prob.value) #结果： 可行解为： [0.45748689] [0.12414101] [0.00955066] [0.12427119] [0.19787315] [0.0866771] 最小值: 0.009432199759803417最小 风险为0.94%
- matlab 二次规划标准型：
m i n 1 2 x T H x + f T x . s . t . { A x ≤ b A e q ⋅ x = b e q l b ≤ x ≤ u b 。 min~\frac{1}{2}x^THx+f^Tx.\\s.t. \left\{\begin{aligned} & Ax\le b\\ &Aeq\cdot x=beq\\ &lb\le x\le ubpython流星雨特效代码简单。 \end{aligned} \right. min 21xTHx+fTx.s.t.⎩ ⎨ ⎧Ax≤bAeq⋅x=beqlb≤x≤ub。
%这里是纯二次型 %注意是2*H [x,value]=quadprog(2*H,zeros(6,1),A,b,Aeq,beq,zeros(6,1))%H为协方差矩阵 %结果： x = 0.4575 0.1241 0.0096 0.1243 0.1979 0.0867 value = 0.0094最小风险为9.4%
- 该二次规划的二次项过多，lingo求解不如matlab和python（矩阵思想）直观，所以这里不提供lingo求解
-
5.4 某糖果厂用原料A,B,C加工成三种不同牌号的糖果甲，乙，丙。已知各种牌号糖果中A,B,C含量，原料成本，各种原料的每月限制用量，三种牌号的糖果的单位加工费以及售价，如表5.7所示。问该厂每月生产者三种牌号糖果各多少千克，才能使其获利最大。试建立这个问题的线性规划问题模型。
原料 甲 乙 丙 原料成本/（元/kg） 每月限制用量/kg A ≥ \ge ≥ 60% ≥ \ge ≥ 30% 2 2000 B 1.5 2500 C ≤ \le ≤ 20% ≤ \le ≤ 50% ≤ \le ≤ 60% 1 1200 加工费/（元/kg） 0.5 0.4 0.3 售价/（元/kg） 3.4 2.85 2.25 -
决策变量： x 1 , x 2 , x 3 ， x 4 , x 5 , x 6 x_1,x_2,x_3，x_4,x_5,x_6 x1,x2,x3，x4,x5,x6 （甲，乙，丙三种糖果在各自每种原料上的质量）
原料 甲 乙 丙 A x1 x2 x3 B x4 x5 x6 C x7 x8 x9 -
利润：售价-加工费-原料成本
-
目标函数 z = 2.9 ( x 1 + x 4 + x 7 ) + 2.45 ( x 2 + x 5 + x 8 ) + 1.95 ( x 3 + x 6 + x 9 ) − 2 ( x 1 + x 2 + x 3 ) − 1.5 ( x 4 + x 5 + x 6 ) − ( x 7 + x 8 + x 9 ) z=2.9(x_1+x_4+x_7)+2.45(x_2+x_5+x_8)+1.95(x_3+x_6+x_9)-2(x_1+x_2+x_3)-1.5(x_4+x_5+x_6)\\-(x_7+x_8+x_9) z=2.9(x1+x4+x7)+2.45(x2+x5+x8)+1.95(x3+x6+x9)−2(x1+x2+x3)−1.5(x4+x5+x6)−(x7+x8+x9)
-
约束条件：
- A，B, C用料限制：
- A： x 1 + x 2 + x 3 ≤ 2000 x_1+x_2+x_3\le2000 x1+x2+x3≤2000
- B: x 4 + x 5 + x 6 ≤ 2500 x_4+x_5+x_6\le2500 x4+x5+x6≤2500
- C: x 7 + x 8 + x 9 ≤ 1200 x_7+x_8+x_9\le1200 x7+x8+x9≤1200
- 甲，乙，丙中A,B,C的含量限制：
- 甲： x 1 x 1 + x 4 + x 7 ≥ 60 % , x 7 x 1 + x 4 + x 7 ≤ 20 % \frac{x_1}{x_1+x_4+x_7}\ge60\%,\frac{x_7}{x_1+x_4+x_7}\le20\% x1+x4+x7x1≥60%,x1+x4+x7x7≤20%
- 乙： x 2 x 2 + x 5 + x 8 ≥ 30 % , x 8 x 2 + x 5 + x 8 ≤ 50 % \frac{x_2}{x_2+x_5+x_8}\ge30\%,\frac{x_8}{x_2+x_5+x_8}\le50\% x2+x5+x8x2≥30%,x2+x5+x8x8≤50%
- 丙： x 9 x 3 + x 6 + x 9 ≤ 60 % \frac{x_9}{x_3+x_6+x_9}\le60\% x3+x6+x9x9≤60%
- A，B, C用料限制：
-
建立线性规划模型：
m a x z = [ 0.9 , 0.45 , − 0.05 , 1.4 , 0.95 , 0.45 , 1.9 , 1.45 , 0.95 ] ⋅ X T s . t . { x 1 + x 2 + x 3 ≤ 2000 x 4 + x 5 + x 6 ≤ 2500 x 7 + x 8 + x 9 ≤ 1200 0.6 ( x 1 + x 4 + x 7 ) ≤ x 1 − 0.2 ( x 1 + x 4 + x 7 ) ≤ − x 7 0.3 ( x 2 + x 5 + x 8 ) ≤ x 2 − 0.5 ( x 2 + x 5 + x 8 ) ≤ − x 8 − 0.6 ( x 3 + x 6 + x 9 ) ≤ x 9 max ~z=[0.9,0.45,-0.05,1.4,0.95,0.45,1.9,1.45,0.95]\cdot X^T\\s.t. \left\{\begin{aligned} &x_1+x_2+x_3\le2000\\ &x_4+x_5+x_6\le2500\\ &x_7+x_8+x_9\le1200\\ &0.6(x_1+x_4+x_7)\le x_1\\ &-0.2(x_1+x_4+x_7)\le -x_7\\ &0.3(x_2+x_5+x_8)\le x_2\\ &-0.5(x_2+x_5+x_8)\le -x_8\\ &-0.6(x_3+x_6+x_9)\le x_9 \end{aligned} \right. max z=[0.9,0.45,−0.05,1.4,0.95,0.45,1.9,1.45,0.95]⋅XTs.t.⎩ ⎨ ⎧x1+x2+x3≤2000x4+x5+x6≤2500x7+x8+x9≤12000.6(x1+x4+x7)≤x1−0.2(x1+x4+x7)≤−x70.3(x2+x5+x8)≤x2−0.5(x2+x5+x8)≤−x8−0.6(x3+x6+x9)≤x9- python
import numpy as np import cvxpy as cp c=np.array([0.9,0.45,-0.05,1.4,0.95,0.45,1.9,1.45,0.95]).reshape(9,1) x=cp.Variable((1,9)) obj=cp.Maximize(x@c) con=[x[0,0]+x[0,1]+[0,2]<=2000,x[0,3]+x[0,4]+x[0,5]<=2500,x[0,6]+x[0,7]+x[0,8]<=1200, 0.6*(x[0,0]+x[0,3]+x[0,6])<=x[0,0],-0.2*(x[0,0]+x[0,3]+x[0,6])<=-x[0,6], 0.3*(x[0,1]+x[0,4]+x[0,7])<=x[0,1],-0.5*(x[0,1]+x[0,4]+x[0,7])<=-x[0,7], -0.6*(x[0,2]+x[0,5]+x[0,8])<=x[0,8],x>=np.zeros((1,9)) ] prob=cp.Problem(obj,con) prob.solve(cp.SCS,verbose=True) print(x.value) print(prob.value) #求解出来的最大利润5446.98650061946- lingo
sets: fac/1..9/:c,x; endsets data: c=0.9,0.45,-0.05,1.4,0.95,0.45,1.9,1.45,0.95; enddata max =@sum(fac(i):c(i)*x(i)); x(1)+x(2)+x(3)<2000; x(4)+x(5)+x(6)<2500; x(7)+x(8)+x(9)<1200; 0.6*(x(1)+x(4)+x(7))<x(1); -0.2*(x(1)+x(4)+x(7))<-x(7); 0.3*(x(2)+x(5)+x(8))<x(2); -0.5*(x(2)+x(5)+x(8))<-x(8); -0.6*(x(3)+x(6)+x(9))<x(9); #结果 X( 1) 580.0000 X( 2) 1420.000 X( 3) 0.000000 X( 4) 386.6667 X( 5) 2113.333 X( 6) 0.000000 X( 7) 0.000000 X( 8) 1200.000 X( 9) 0.000000可以看到lingo的代码非常直观
最大利润：5450
- matlab 在这里构造系数矩阵A比较繁琐，不考虑matlab实现（不够直观）
-
-
5.5 试求解多目标线性规划问题：
m a x z 1 = 100 x 1 + 90 x 2 + 80 x 3 + 70 x 4 m i n z 2 = 3 x 2 + 2 x 4 { x 1 + x 2 ≥ 30 x 3 + x 4 ≥ 30 3 x 1 + 2 x 3 ≤ 120 3 x 2 + 2 x 4 ≤ 48 x i ≥ 0 , i = 1 , 2 , 3 , 4 \begin{aligned} &max~~z_1=100x_1+90x_2+80x_3+70x_4\\ &min~~z_2=3x_2+2x_4\\ &\left\{\begin{aligned} &x_1+x_2\ge30\\ &x_3+x_4\ge 30\\ &3x_1+2x_3\le120\\ &3x_2+2x_4\le48\\ &x_i\ge0,~i=1,2,3,4 \end{aligned}\right. \end{aligned} max z1=100x1+90x2+80x3+70x4min z2=3x2+2x4⎩ ⎨ ⎧x1+x2≥30x3+x4≥303x1+2x3≤1203x2+2x4≤48xi≥0, i=1,2,3,4- 求解该多目标问题，只需要构造 f = m i n ( z 2 − z 1 ) f=min(z_2-z_1) f=min(z2−z1)即可
- m i n f = [ − 100 , − 87 , − 80 , − 68 ] ⋅ X T min ~f=[-100,-87,-80,-68]\cdot X^T min f=[−100,−87,−80,−68]⋅XT
- matlab
c=[-100,-87,-80,-68]; A=[-1,-1,0,0;0,0,-1,-1;3,0,2,0;0,3,0,2]; b=[-30,-30,120,48]'; [x,fval]=linprog(c,A,b,[],[],zeros(4,1)); %结果 fval = -5.9120e+03 >> x x = 14.0000 16.0000 39.0000 0- python
import cvxpy as cp import numpy as np c=np.array([-100,-87,-80,-68]).reshape((1,4)) x=cp.Variable((4,1)) b=np.array([-30,-30,120,48]).reshape(4,1) A=np.array([[-1,-1,0,0],[0,0,-1,-1],[3,0,2,0],[0,3,0,2]]) obj=cp.Minimize(c@x) con=[A@x<=b,x>=np.zeros((4,1))] pro=cp.Problem(obj,con) pro.solve() print(pro.value) print(x.value) #结果： -5911.999999475287 [[1.40000001e+01] [1.59999999e+01] [3.89999998e+01] [1.26453538e-07]]- lingo 在这里反而还需要将不等式一一输入，相较于矩阵的思想变得繁琐（不做赘述）
综上可以发现求解线性规划问题没有最好的方法，那种方法更直观就使用哪种方法
2901
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
