devbox
我家小花儿
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

数学建模——线性规划类_python数学建模线性规划代码

作者：我家小花儿 | 2024-08-06 17:53:41

赞

踩

python数学建模线性规划代码

一.线性规划

1.算法

【1】通用代码

[x,y]=linprog(c,A,b,Aeq,beq,lb,ub)

例如：

max需要加负号变成min、>=需要加负号变成<=

matlab

（1）基于求解器

  1. %    LP1_1_1.m
  2.  
  3. f=[-2;-3;5];
  4. a=[-2,5,-1;1,3,1];b=[-10;12];
  5. aeq=[1,1,1];beq=7;
  6. [x,y]=linprog(f,a,b,aeq,beq,zeros(3,1));
  7. x,y=-y    %最大值
  8.  
  9.  
  10. x =
  11.     6.4286
  12.     0.5714
  13.          0
  14. y =
  15.    14.5714
  16. % 说明：在x1 x2 x3 = 6.4286 0.5714 0 的情况下，maxz = 14.5714

（2）基于问题

con中根据符号分类

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

python

  1. #    线性规划模型.py
  2.  
  3. from scipy.optimize import linprog
  4. import numpy as np
  5.  
  6. c=np.array([-2,-3,5])
  7. a=[[-2,5,-1],[1,3,1]]
  8. b=[[-10],[12]]
  9. aeq=[[1,1,1]]
  10. beq=[7]
  11. lb=[0,0,0]
  12. ub=[None]*len(c)
  13. bound=tuple(zip(lb,ub))
  14. res=linprog(c,a,b,aeq,beq,bound,method='simplex',options={"disp":True})
  15. print("最优解：",res.x)
  16. print("目标函数最小值：",res.fun)
  17. print("目标函数最大值：",-res.fun)
  18. print(res)
  19.  
  20. 最优解： [6.42857143 0.57142857 0.        ]
  21. 目标函数最小值： -14.571428571428571
  22. 目标函数最大值： 14.571428571428571
  23.      con: array([0.])
  24.      fun: -14.571428571428571
  25.  message: 'Optimization terminated successfully.'
  26.      nit: 3
  27.    slack: array([0.        , 3.85714286])
  28.   status: 0
  29.  success: True
  30.        x: array([6.42857143, 0.57142857, 0.        ])

【2】可以转化为线性规划

（1）绝对值

  1. %    LP_abs.m
  2.  
  3. clc,clear
  4. c=[1,2,3,4]';
  5. a=[1,-1,-1,1;1,-1,1,-3;1,-1,-2,3];
  6. b=[-2,-1,-1/2]';
  7. prob=optimproblem;
  8. u=optimvar('u',4,'LowerBound',0);
  9. v=optimvar('v',4,'LowerBound',0);
  10. prob.Objective=sum(c'*(u+v));
  11. prob.Constraints.con=a*(u-v)<=b;
  12. [sol,fval,flag,out]=solve(prob)
  13. x=sol.u-sol.v,minz=fval
  14.  
  15.  
  16. x =
  17.     -2
  18.      0
  19.      0
  20.      0
  21. minz =
  22.     2.0000

（2）min(max(q*x))

（见风投案例模型二）

2.风投案例

【0】题目描述

【1】模型一

模型一：设定风险度的最大接受值，在不太冒险的情况下选择最大收益。

（1）matlab

  1. %    LP_fengtou_1.m
  2.  
  3. clc,clear
  4. t=0;hold on
  5. while t<0.05
  6.     c=[-0.05,-0.27,-0.19,-0.185,-0.185];
  7.     A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];
  8.     b=t*ones(4,1);
  9.     Aeq=[1,1.01,1.02,1.045,1.065];beq=1;
  10.     lb=zeros(5,1);
  11.     [x,Q]=linprog(c,A,b,Aeq,beq,lb);
  12.     Q=-Q;
  13.     plot(t,Q,'*k');
  14.     t=t+0.001;
  15. end
  16. xlabel('t'),ylabel('Q')

选择曲线的转折点作为最优投资组合。

输出转折点见模型一python代码：

鼠标在图上确定大致坐标，在循环中加入条件输出。

（2）python

  1. #    风投案例.py
  2.  
  3. from scipy.optimize import linprog
  4. import numpy as np
  5. from numpy import ones,diag,c_,zeros
  6. import matplotlib.pyplot as plt
  7.  
  8. plt.rc('font',size=16)
  9. c=np.array([-0.05,-0.27,-0.19,-0.185,-0.185])
  10. A=c_[zeros(4),diag([0.025,0.015,0.055,0.026])]
  11. Aeq=[[1,1.01,1.02,1.045,1.065]]
  12. beq=[1]
  13. a=0
  14. aa=[]
  15. ss=[]
  16. while a<0.05:
  17.     b=ones(4)*a
  18.     res = linprog(c, A, b, Aeq, beq)
  19.     x=res.x
  20.     Q=-res.fun
  21.     aa.append(a)
  22.     ss.append(Q)
  23. #输出转折点
  24.     if a==0.006:
  25.         print("x=", x)
  26.         print("Q=", Q)
  27.     a=a+0.001
  28. plt.plot(aa,ss,'r*')
  29. plt.xlabel('$a$')
  30. plt.ylabel('$Q$',rotation=90)
  31. plt.savefig('模型一.png',dpi=500)
  32. plt.show()
  33.  
  34. x= [0.         0.24       0.4        0.10909091 0.22122066]
  35. Q= 0.20190763977806234（万元）

【2】模型二

模型二的目标函数：使各项目中风险最大的一项风险最小。

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import cvxpy as cp
  4.  
  5. plt.rc('font',size=16)
  6. x=cp.Variable(6,pos=True)
  7. obj=cp.Minimize(x[5])
  8. a1=np.array([0.025,0.015,0.055,0.026])
  9. a2=np.array([0.05,0.27,0.19,0.185,0.185])
  10. a3=np.array([1,1.01,1.02,1.045,1.065])
  11. k=0.05
  12. kk=[]
  13. ss=[]
  14. while k<0.27:
  15.     con=[cp.multiply(a1,x[1:5])-x[5]<=0,
  16.          0.05*x[0]+0.27*x[1]+0.19*x[2]+0.185*x[3]+0.185*x[4]>=k,
  17.          x[0]+1.01*x[1]+1.02*x[2]+1.045*x[3]+1.065*x[4]==1]
  18.     prob=cp.Problem(obj,con)
  19.     prob.solve(solver='CVXOPT')
  20. #输出转折点
  21.     if abs(k-0.21)<0.0005:
  22.         print("x=",x.value)
  23.         print("风险=",prob.value) 
  24.     kk.append(k)
  25.     ss.append(prob.value)
  26.     k=k+0.005
  27. plt.plot(kk,ss,'r*')
  28. plt.xlabel('$k$')
  29. plt.ylabel('$R$',rotation=90)
  30. plt.savefig('模型二.png',dpi=500)
  31. plt.show()
  32.  
  33. x= [2.09147087e-08 3.08874349e-01 5.14790372e-01 1.40396702e-01
  34.  1.52452147e-02 7.72185738e-03]
  35. 风险= 0.00772185738406626

【3】模型三

模型三：风险与收益的加权线性组合。

（1）matlab

  1. clc,clear
  2. M=10000;prob=optimproblem;
  3. x=optimvar('x',6,1,'LowerBound',0);
  4. r=[0.05,0.28,0.21,0.23,0.25];
  5. p=[0,0.01,0.02,0.045,0.065];
  6. q=[0,0.025,0.015,0.055,0.026]';
  7. %w=0:0.1:1;
  8. w=[0.766,0.767,0.810,0.811,0.824,0.825,0.962,0.963,1.0];
  9. V=[];%风险
  10. Q=[];%收益
  11. X=[];%最优解
  12. prob.Constraints.con1=(1+p)*x(1:end-1)==M;
  13. prob.Constraints.con2=q(2:end).*x(2:end-1)<=x(end);
  14. for i=1:length(w)
  15.     prob.Objective=w(i)*x(end)-(1-w(i))*(r-p)*x(1:end-1);
  16.     [sol,fval,flag,out]=solve(prob);
  17.     xx=sol.x;
  18.     V=[V,max(q.*xx(1:end-1))];
  19.     Q=[Q,(r-p)*xx(1:end-1)];
  20.     X=[X;xx'];
  21.     plot(V,Q,'*-');grid on
  22.     xlabel('风险');ylabel('收益')
  23.     X
  24. end
  25. V,Q,format
  26.  
  27. w=（权重）
  28.     [0.766,    0.767,    0.810,    0.811,    0.824,    0.825,    0.962,    0.963,    1.0];
  29. V =（风险）
  30.   247.5248   92.2509   92.2509   78.4929   78.4929   59.3960   59.3960         0         0
  31. Q =（收益）
  32.    1.0e+03 *
  33.     2.6733    2.1648    2.1648    2.1060    2.1060    2.0162    2.0162    0.5000    0.5000
  34. X =
  35.    1.0e+04 *
  36.  
  37.          0    0.9901         0         0         0    0.0248
  38.          0    0.3690    0.6150         0         0    0.0092
  39.          0    0.3690    0.6150         0         0    0.0092
  40.          0    0.3140    0.5233    0.1427         0    0.0078
  41.          0    0.3140    0.5233    0.1427         0    0.0078
  42.          0    0.2376    0.3960    0.1080    0.2284    0.0059（转折点处的X值）
  43.          0    0.2376    0.3960    0.1080    0.2284    0.0059
  44.     1.0000         0         0         0         0         0
  45.     1.0000         0         0         0         0         0

（2）python

  1. from scipy.optimize import linprog
  2. import matplotlib.pyplot as plt
  3.  
  4. plt.rc('font',size=16)
  5. A=[[0,0.025,0,0,0,-1],[0,0,0.15,0,0,-1],[0,0,0,0.55,0,-1],[0,0,0,0,0.026,-1]]
  6. b=[0,0,0,0]
  7. Aeq=[[1,1.01,1.02,1.045,1.065,0]]
  8. beq=[1]
  9. bound=((0,None),(0,None),(0,None),(0,None),(0,None),(0,None))
  10. s=0;ss=[];aa=[]
  11. while s<=1:
  12.     c=[-(1-s)*0.05,-(1-s)*0.27,-(1-s)*0.19,-(1-s)*0.185,-(1-s)*0.185,s]
  13.     res=linprog(c,A,b,Aeq,beq)
  14.     ss.append(s);aa.append(-res.fun)
  15.     s=s+0.01
  16. plt.plot(ss,aa,'r.')
  17. plt.show()
  1. import matplotlib.pyplot as plt
  2. import cvxpy as cp
  3. import numpy as np
  4.  
  5. plt.rc('font',size=16)
  6. x=cp.Variable(6)
  7. a1=np.array([0.025,0.015,0.055,0.026])
  8. a2=np.array([0.05,0.27,0.19,0.185,0.185])
  9. a3=np.array([1,1.01,1.02,1.045,1.065])
  10. con=[cp.multiply(a1,x[1:5])-x[5]<=0,a3@x[:-1]==1,x>=0]
  11. s=0;ss=[];aa=[]
  12. while s<=1:
  13.     obj=cp.Minimize(s*x[-1]-(1-s)*a2@x[:-1])
  14.     prob=cp.Problem(obj,con)
  15.     prob.solve(solver='GLPK_MI')
  16.     ss.append(s);aa.append(-prob.value)
  17.     s=s+0.01
  18. plt.plot(ss,aa,'r.')
  19. plt.show()

二.整数规划

0.分类

1.通用代码

  1. %    int_LP_1.m
  2.  
  3. clc,clear,prob=optimproblem;
  4. x=optimvar('x',6,'Type','integer','LowerBound',0);
  5. prob.Objective=sum(x);
  6. % 创建一个由空优化约束组成的 6×1 数组，使用 constr 初始化用于创建约束表达式的循环
  7. coon=optimconstr(6);
  8. a=[35,40,50,45,55,30];
  9. con(1)=x(1)+x(6)>=35;
  10. for i=1:5
  11.     con(i+1)=x(i)+x(i+1)>=a(i+1);
  12. end
  13. prob.Constraints.con=con;
  14. [sol,fval,flag]=solve(prob)
  15. x=sol.x
  16.  
  17. x =
  18.     35
  19.      5
  20.     45
  21.      0
  22.     55
  23.      0

intlinprog函数

x = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)：函数相比linprog多了一个参数intcon，用于标定整数变量的位置：x1、x2为整数，即intcon = [1 2]

  1. %    int_LP_intlinprog.m
  2.  
  3. c = [-40;-90];
  4. A = [9,7;7,20];
  5. b = [56;70];
  6. [x,fval] = intlinprog(c,[1 2],A,b,[],[],zeros(2,1));
  7. x,maxz=-fval
  8.  
  9. x =
  10.     4.0000
  11.     2.0000
  12. maxz =
  13.    340

2.分支定界法（剪枝）

尝试一维单调

  1. %    int_branchbunch.m
  2. clear global;
  3. clear;
  4. clc;
  5.  
  6. global result; % 存储所有整数解
  7. global lowerBound; % 下界
  8. global upperBound; % 上界
  9. global count; % 用以判断是否为第一次分支
  10.  
  11. count = 1;
  12.  
  13. f = [-40, -90];
  14. A = [8, 7;7, 20;];  
  15. b = [56; 70];
  16. Aeq = [];
  17. beq = [];
  18. lbnd = [0; 0];
  19. ubnd = [inf; inf];
  20.  
  21. % f = [-3 2 -5];
  22. % A = [1 2 -1;1 4 1;1 1 0;0 4 1];
  23. % b = [2;4;3;6];
  24. % Aeq = [];
  25. % beq = [];
  26. % lbnd=[0 0 0];
  27. % ubnd=[1 1 1];
  28.  
  29. BinTree = createBinTreeNode({f, A, b, Aeq, beq, lbnd, ubnd});
  30. if ~isempty(result)
  31.     [fval,flag]=min(result(:,end)); % result中每一行对应一个整数解及对应的函数值
  32.     Result=result(flag,:);
  33.     disp('该整数规划问题的最优解为：');
  34.     disp(Result(1,1:end-1));
  35.     disp('该整数规划问题的最优值为：');
  36.     disp(Result(1,end));
  37. else
  38.     disp('该整数规划问题无可行解');
  39. end
  40.  
  41. function BinTree = createBinTreeNode(binTreeNode)
  42.  
  43. global result;
  44. global lowerBound;
  45. global upperBound;
  46. global count;
  47.  
  48. if isempty(binTreeNode)
  49.     return;
  50. else
  51.     BinTree{1} = binTreeNode;
  52.     BinTree{2} = [];
  53.     BinTree{3} = [];
  54.     
  55.     [x, fval, exitflag] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  56.         binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, binTreeNode{7});
  57.     if count == 1
  58. %         upperBound = 0; % 初始下界为空
  59.         lowerBound = fval;
  60.         count = 2;
  61.     end
  62.     
  63.     if ~IsInRange(fval)
  64.         return;
  65.     end
  66.  
  67.     if exitflag == 1
  68.         flag = [];
  69.         % 寻找非整数解分量
  70.         for i = 1 : length(x)
  71.             if round(x(i)) ~= x(i)
  72.                 flag = i;
  73.                 break;
  74.             end
  75.         end
  76.         % 分支
  77.         if ~isempty(flag)
  78.             lowerBound = max([lowerBound; fval]);
  79.             OutputLowerAndUpperBounds();
  80.             lbnd = binTreeNode{6};
  81.             ubnd = binTreeNode{7};
  82.             lbnd(flag, 1) = ceil(x(flag, 1)); % 朝正无穷四舍五入
  83.             ubnd(flag, 1) = floor(x(flag, 1));
  84.             
  85.             % 进行比较，优先选择目标函数较大的进行分支
  86.             [~, fval1] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  87.         binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd);
  88.             [~, fval2] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  89.         binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7});
  90.             if fval1 < fval2                
  91.                 % 创建左子树          
  92.                 BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  93.             binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
  94.  
  95.                 % 创建右子树
  96.                 BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  97.             binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
  98.             else
  99.                 % 创建右子树
  100.                 BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  101.             binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
  102.                 % 创建左子树          
  103.                 BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
  104.             binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
  105.             end
  106.         else
  107.             upperBound = min([upperBound; fval]);
  108.             OutputLowerAndUpperBounds();
  109.             result = [result; [x', fval]];
  110.             return;
  111.         end
  112.     else
  113.         % 剪枝
  114.         return;
  115.     end  
  116. end
  117. end
  118.  
  119. function y = IsInRange(fval)
  120.     global lowerBound;
  121.     global upperBound;
  122.  
  123.     if isempty(upperBound) & fval >= lowerBound
  124.         y = 1;
  125.     else if  (fval >= lowerBound & fval <= upperBound)
  126.         y = 1;
  127.     else
  128.         y = 0;
  129.     end
  130. end
  131. end
  132.  
  133. function y = OutputLowerAndUpperBounds()
  134.  
  135. global lowerBound;
  136. global upperBound;
  137.  
  138. disp("此时下界、上界分别为");
  139. disp(lowerBound);
  140. if isempty(upperBound)
  141.     disp('  无上界')
  142. else
  143.     disp(upperBound);
  144. end
  145.  
  146. end
  147.  
  148. 该整数规划问题的最优解为：
  149.      4     2
  150. 该整数规划问题的最优值为：
  151.   -340

3.割平面法

割去小数部分，尝试二维单调

  1. function  [intx,intf] = DividePlane(A,c,b,baseVector)
  2. %功能：用割平面法求解整数规划
  3. %调用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
  4. %其中, A：约束矩阵；
  5. %      c：目标函数系数向量；
  6. %      b：约束右端向量；
  7. %      baseVector：初始基向量；
  8. %      intx：目标函数取最值时的自变量值；
  9. %      intf：目标函数的最值；
  10. sz = size(A);
  11. nVia = sz(2);%获取有多少决策变量
  12. n = sz(1);%获取有多少约束条件
  13. xx = 1:nVia;
  14.  
  15. if length(baseVector) ~= n
  16.     disp('基变量的个数要与约束矩阵的行数相等！');
  17.     mx = NaN;
  18.     mf = NaN;
  19.     return;
  20. end
  21.  
  22. M = 0;
  23. sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
  24. xb = b;
  25.  
  26. %首先用单纯形法求出最优解
  27. while 1   
  28.     [maxs,ind] = max(sigma);
  29. %--------------------用单纯形法求最优解--------------------------------------
  30.     if maxs <= 0   %当检验数均小于0时，求得最优解。      
  31.         vr = find(c~=0 ,1,'last');
  32.         for l=1:vr
  33.             ele = find(baseVector == l,1);
  34.             if(isempty(ele))
  35.                 mx(l) = 0;
  36.             else
  37.                 mx(l)=xb(ele);
  38.             end
  39.         end
  40.         if max(abs(round(mx) - mx))<1.0e-7  %判断最优解是否为整数解，如果是整数解。
  41.             intx = mx;
  42.             intf = mx*c;
  43.             return;
  44.         else  %如果最优解不是整数解时，构建切割方程
  45.             sz = size(A);
  46.             sr = sz(1);
  47.             sc = sz(2);
  48.             [max_x, index_x] = max(abs(round(mx) - mx));
  49.             [isB, num] = find(index_x == baseVector);
  50.             fi = xb(num) - floor(xb(num));
  51.             for i=1:(index_x-1)
  52.                 Atmp(1,i) = A(num,i) - floor(A(num,i));
  53.             end
  54.             for i=(index_x+1):sc
  55.                 Atmp(1,i) = A(num,i) - floor(A(num,i));
  56.             end
  57.             
  58.             Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
  59.             A = [A zeros(sr,1);-Atmp(1,:) 1];
  60.             xb = [xb;-fi];
  61.             baseVector = [baseVector sc+1];
  62.             sigma = [sigma 0];
  63.          
  64.             %-------------------对偶单纯形法的迭代过程----------------------
  65.             while 1
  66.                 %----------------------------------------------------------
  67.                 if xb >= 0    %判断如果右端向量均大于0，求得最优解
  68.                     if max(abs(round(xb) - xb))<1.0e-7   
  69. %如果用对偶单纯形法求得了整数解，则返回最优整数解
  70.                         vr = find(c~=0 ,1,'last');
  71.                         for l=1:vr
  72.                             ele = find(baseVector == l,1);
  73.                             if(isempty(ele))
  74.                                 mx_1(l) = 0;
  75.                             else
  76.                                 mx_1(l)=xb(ele);
  77.                             end
  78.                         end
  79.                         intx = mx_1;
  80.                         intf = mx_1*c;
  81.                         return;
  82.                     else   %如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程
  83.                         sz = size(A);
  84.                         sr = sz(1);
  85.                         sc = sz(2);
  86.                         [max_x, index_x] = max(abs(round(mx_1) - mx_1));
  87.                         [isB, num] = find(index_x == baseVector);
  88.                         fi = xb(num) - floor(xb(num));
  89.                         for i=1:(index_x-1)
  90.                             Atmp(1,i) = A(num,i) - floor(A(num,i));
  91.                         end
  92.                         for i=(index_x+1):sc
  93.                             Atmp(1,i) = A(num,i) - floor(A(num,i));
  94.                         end
  95.                         Atmp(1,index_x) = 0;  %下一次对偶单纯形迭代的初始表格
  96.                         A = [A zeros(sr,1);-Atmp(1,:) 1];
  97.                         xb = [xb;-fi];
  98.                         baseVector = [baseVector sc+1];
  99.                         sigma = [sigma 0];
  100.                         continue;
  101.                     end
  102.                 else   %如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
  103.                     minb_1 = inf;
  104.                     chagB_1 = inf;
  105.                     sA = size(A);
  106.                     [br,idb] = min(xb);
  107.                     for j=1:sA(2)
  108.                         if A(idb,j)<0
  109.                             bm = sigma(j)/A(idb,j);
  110.                             if bm<minb_1
  111.                                 minb_1 = bm;
  112.                                 chagB_1 = j;
  113.                             end
  114.                         end
  115.                     end
  116.                     sigma = sigma -A(idb,:)*minb_1;
  117.                     xb(idb) = xb(idb)/A(idb,chagB_1);
  118.                     A(idb,:) = A(idb,:)/A(idb,chagB_1);
  119.                     for i =1:sA(1)
  120.                         if i ~= idb
  121.                             xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
  122.                             A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
  123.                         end
  124.                     end
  125.                     baseVector(idb) = chagB_1;
  126.                 end
  127.               %------------------------------------------------------------
  128.             end 
  129.             %--------------------对偶单纯形法的迭代过程---------------------    
  130.         end     
  131.     else     %如果检验数有不小于0的，则进行单纯形算法的迭代过程
  132.         minb = inf;
  133.         chagB = inf;
  134.         for j=1:n
  135.             if A(j,ind)>0
  136.                 bz = xb(j)/A(j,ind);
  137.                 if bz<minb
  138.                     minb = bz;
  139.                     chagB = j;
  140.                 end
  141.             end
  142.         end
  143.         sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
  144.         xb(chagB) = xb(chagB)/A(chagB,ind);
  145.         A(chagB,:) = A(chagB,:)/A(chagB,ind);
  146.         for i =1:n
  147.             if i ~= chagB
  148.                 xb(i) = xb(i)-A(i,ind)*xb(chagB);
  149.                 A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
  150.             end
  151.         end
  152.         baseVector(chagB) = ind;
  153.     end
  154.     M = M + 1;
  155.     if (M == 1000000)
  156.         disp('找不到最优解！');
  157.         mx = NaN; 
  158.         minf = NaN;
  159.         return;
  160.     end
  161. end
  1. %    int_LP_divideplane.m
  2.  
  3. c = [-1;-1]; % 不要加松弛变量
  4. A = [-1 1 1 0;3 1 0 1]; % 加上松弛变量
  5. b = [1;4];
  6. [x,fval] = DividePlane(A,c,b,[3 4]); % 松弛变量 3 4
  7. x,maxz=fval
  8.  
  9. x =
  10.     1.0000    1.0000
  11. maxz =
  12.     -2

4.匈牙利算法

【1】

【2】

  1. %    int_XiongYaLi.m
  2.  
  3. c = [6 7 11 2;4 5 9 8;3 1 10 4;5 9 8 2];
  4. a = zeros(8,16);
  5. for i = 1:4
  6. a(i,(i-1)*4+1:4*i) = 1;
  7. a(4+i,i:4:16) = 1;
  8. end
  9. b = ones(8,1);
  10. [x,y] = linprog(c,[],[],a,b,zeros(16,1),ones(16,1));
  11. X = reshape(x,4,4)
  12. opt = y
  13.  
  14. % X =
  15. % 
  16. %      0     0     0     1
  17. %      1     0     0     0
  18. %      0     1     0     0
  19. %      0     0     1     0
  20. % 
  21. % 
  22. % opt =
  23. % 
  24. %     15
  25.  
  26. A=[6 7 11 2;
  27.     4 5 9 8;
  28.     3 1 10 4;
  29.     5 9 8 2];
  30. % A = [3 8 2 10 3;8 7 2 9 7;6 4 2 7 5;8 4 2 3 5;9 10 6 9 10];
  31. B=Hungary(A);
  32. [~,b]=linear_assignment(A,B)
  33.  
  34. % b =
  35. %        4              1              8              2 
  36.  
  37. function res=Hungary(N)
  38. %输入的矩阵应N*N的
  39. [a,~]=size(N);
  40. %第一步每一行减去当前行最小值
  41. for ii = 1:a
  42.     N(ii,:)= N(ii,:)-min( N(ii,:));
  43. end
  44. %第二步每一列减去当前列最小值
  45. for ii = 1:a
  46.     N(:,ii)=  N(:,ii)-min( N(:,ii));
  47. end
  48. num=0;
  49. while num~=a
  50.     [num,N_min,del_hang,del_lie]=line_count(N);
  51.     if num ~=a
  52.         for ii=1:a
  53.             if del_hang(ii)~=ii
  54.                 N(ii,:) =  N(ii,:)-N_min;
  55.             end
  56.             if del_lie(ii)==ii
  57.             N(:,ii) =  N(:,ii)+N_min;
  58.             end
  59.         end
  60.     else
  61.         res=N;
  62.     end
  63. end
  64.  
  65. function [num,M_min,del_hang,del_lie]=line_count(M)
  66. [a,~]=size(M);
  67. num=0;
  68. h=0;
  69. del_hang=zeros(a,1);
  70. del_lie=zeros(a,1);
  71. for ii=1:a
  72.     del=ii-h;
  73.     [~,b]=size(find(M(del,:)==0));
  74.     if   b>= 2
  75.         M(del,:)=[];
  76.         h=h+1;
  77.         del_hang(ii)=ii;    %得到被覆盖的行数
  78.         num=num+1;
  79.     end
  80. end
  81. l=0;
  82. for ii=1:a
  83.     del=ii-l;
  84.     [b,~]=size(find(M(:,del)==0));
  85.     if  b >=1
  86.         M(:,del)=[];
  87.         l=l+1;
  88.         del_lie(ii)=ii;    %得到被覆盖的列数
  89.         num=num+1;
  90.     end
  91. end
  92. M_min=min(min(M));
  93. end
  94. end
  95.  
  96. function [place,res]=linear_assignment(M,N)
  97. %N是n维矩阵,N是经过Hungary处理的
  98. %M是未处理前的
  99. [a,~]=size(N);
  100. x=0;
  101. place=zeros(1,a);
  102. res=zeros(1,a);
  103. judge=zeros(1,a);
  104. while find(N==0)
  105.     for ii=1:a
  106.     judge(ii)=length(find(N(ii,:)==0));
  107.     end
  108.     judge(find(judge==0))=[];
  109.     if min(judge)==1
  110.      for ii=1:a
  111.         if length(find(N(ii,:)==0))==1     %先选出行中只有1个0
  112.             x=x+1;
  113.             place(x)=ii+(find(N(ii,:)==0)-1)*a; %得到矩阵中的位置
  114.             h=find(N(ii,:)==0);
  115.             N(ii,:)=1./zeros(1,a);
  116.             N(:,h)=1./zeros(a,1);
  117.         end
  118.      end
  119.     end
  120.     
  121.     for ii=1:a
  122.     judge(ii)=length(find(N(ii,:)==0));
  123.     end
  124.     judge(find(judge==0))=[];
  125.     
  126.     if min(judge)==2
  127.        x=x+1;
  128.     q=find(N==0);
  129.     place(x)=q(1);
  130.     N(mod(q(1),a),:)=1./zeros(1,a);
  131.     N(:,fix(q(1)/a)+1)=1./zeros(a,1);  
  132.     end
  133. end
  134. [place,~]=sort(place);
  135. for ii=1:length(place)
  136.     res(ii)=M(place(ii));
  137. end
  138. end

5.0-1规划

【1】intlinprog函数

0-1规划问题也可以看做区间在[0,1]的整数规划，下面利用intlinprog函数进行计算

  1. %    intlinprog_01.m
  2.  
  3. c = [-3 2 -5]';
  4. A = [1 2 -1;1 4 1;1 1 0;0 4 1];
  5. b = [2;4;3;6];
  6. [x,fval] = intlinprog(c,[1 2 3],A,b,[],[],zeros(3,1),ones(3,1));
  7. x,maxz = -fval
  8.  
  9. x =
  10.      1
  11.      0
  12.      1
  13. maxz =
  14.      8

【2】指派问题

（1）纯0-1整数规划
  1. c=np.array([[4,8,7,15,12],
  2.             [7,9,17,14,10],
  3.             [6,9,12,8,7],
  4.             [6,7,14,6,10],
  5.             [6,9,12,10,6]])
  6. x=cp.Variable((5,5),integer=True)
  7. obj=cp.Minimize(cp.sum(cp.multiply(c,x)))
  8. con=[0<=x,x<=1,
  9.      cp.sum(x,axis=0,keepdims=True)==1,
  10.      cp.sum(x,axis=1,keepdims=True)==1]
  11. prob=cp.Problem(obj,con)
  12. prob.solve(solver='GLPK_MI')
  13. print("最优值为：",prob.value)
  14. print("最优解为：\n",x.value)
  15.  
  16. 最优值为： 34.0
  17. 最优解为：
  18.  [[0. 0. 1. 0. 0.]
  19.  [0. 1. 0. 0. 0.]
  20.  [1. 0. 0. 0. 0.]
  21.  [0. 0. 0. 1. 0.]
  22.  [0. 0. 0. 0. 1.]]
（2）广义指派
  1. 5*x(1)+4*x(2)<=24+y(1)M;        %M是充分大的数
  2. 7*x(1)+3*x(2)<=45+y(2)M;
  3. y(1)+y(2)=1;

i=1,2,3为三种生产方式

x(i)为产量；c(i)为每件产品的成本；k(i)为每种方式的固定成本

  1. min z=(k(1)*y(1)+c(1)*x(1))+(k(2)*y(2)+c(2)*x(2))+(k(3)*y(3)+c(3)*x(3));
  2. y(i)*m<=x(i)<=y(i)*M;   %m充分小，M充分大
  3. y(i)=0 or 1;

【3】旅行商问题(TSP)

（司P31）比赛项目排序问题中，由于开始项目和结束项目没有连接，可引入虚拟项目15，与各项目的权重为0

三.非线性规划

    • 算法

【1】二次规划

目标函数二次，约束条件线性

H为实对称矩阵

[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)
  1. %    DLP_1.m
  2.  
  3. H=[4,-4;-4,8];
  4. f=[-6;3];
  5. a=[1,1;4,1];
  6. b=[3;9];
  7. [x,fval]=quadprog(H,f,a,b,[],[],zeros(2,1))
  8.  
  9. x =
  10.     2.0854
  11.     0.6585
  12. fval =
  13.    -5.5976
  1. %    DLP_2.m
  2.  
  3. x0 = [1;1];
  4. A=[1,1;4,1];
  5. b=[3;9];
  6. Aeq = [];
  7. beq = [];
  8. lb = [0;0];
  9. [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb)
  10. function f = fun(x)
  11.     f = 2*x(1)-4*x(1)*x(2)+4*x(2)^2-6*x(1)-3*x(2)
  12. end
  13.  
  14. x =
  15.     2.0000
  16.     1.0000
  17. fval =
  18.   -15.0000

【2】无约束非线性规划

无约束优化问题有局部最优解的充分条件：梯度=0；Hesse矩阵正定

  1. %    NLP_nocon.m
  2.  
  3. clc,clear
  4. f = @(x) x(1)^3-x(2)^3+3*x(1)^2+3*x(2)^2-9*x(1);
  5. g = @(x) -f(x);
  6. [m1,n1] = fminunc(f,[0,0])%求极小值
  7. [m2,n2] = fminsearch(g,[0,0]);%求极大值
  8. m2,n2=-n2
  9.  
  10. m1 =
  11.     1.0000   -0.0000
  12. n1 =
  13.     -5
  14. m2 =
  15.    -3.0000    2.0000
  16. n2 =
  17.    31.0000

【3】有约束非线性规划

  1. %    NLP_con_1.m
  2. %[x,fval]=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub,'fun2')
  3.  
  4. [x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')
  5.  
  6. function f = fun1(x);
  7. f=sum(x.^2)+8;
  8. end
  9.  
  10. function [g,h] = fun2(x)
  11. g=[-x(1)^2+x(2)-x(3)^2
  12.     x(1)+x(2)^2+x(3)^3-20];
  13. h=[-x(1)-x(2)^2+2
  14.     x(2)+2*x(3)^2-3];
  15. end
  16.  
  17. x =
  18.     0.5522
  19.     1.2033
  20.     0.9478
  21. y =
  22.    10.6511
  1. %    NLP_con_2.m
  2.  
  3. clc,clear
  4. x = optimvar('x',3,'LowerBound',0);
  5. prob = optimproblem('Objective',sum(x.^2)+8);
  6. con1 = [-x(1)^2+x(2)-x(3)^2 <= 0
  7.     x(1)+x(2)^2+x(3)^3 <= 20];
  8. con2 = [-x(1)-x(2)^2+2 == 0
  9.     x(2)+2*x(3)^2 == 3];
  10. prob.Constraints.con1 = con1;
  11. prob.Constraints.con2 = con2;
  12. x0.x = rand(3,1);
  13. [s,f,flag,out] = solve(prob,x0);
  14. s.x,f

（1）只有等式约束

拉格朗日法

（2）一般形式

罚函数法

  1. 不等式约束g(x,i)<=0    等式约束h(x,j)=0
  2. 不等式 g(x,i)<=0  转化为  等式 max(0,g(x,i))=0
  3. 罚函数：T(x.M)=f(x)+M*sum(max(0，g(x,i)))+M*sum(h(x,j)^2)

【4】凸规划（f''>=0）

（1）判定：Hesse处处半正定

  1. %    NLP_TU.m
  2.  
  3. clc,clear
  4. prob=optimproblem;
  5. x=optimvar('x',2,'LowerBound',0);
  6. prob.Objective=sum(x.^2)-4*x(1)+4;
  7. con=[-x(1)+x(2)-2<=0
  8.     x(1)^2-x(2)+1<=0];
  9. prob.Constraints.con=con;
  10. x0.x=rand(2,1)
  11. [sol,fval,flag,out]=solve(prob,x0),sol.x

（2）K-T条件

只要是最优点，一定满足K-T条件

若为凸规划，则充要条件

【5】蒙特卡洛法

  1. %    MTKL.m
  2.  
  3. clc,clear
  4. rng(0)          %rng('shuffle')
  5. p0=0;n=10^6;tic
  6. for i=1:n
  7.     x=randi([0,99],1,5);
  8.     [f,g]=mengte(x);
  9.     if all(g<=0)
  10.         if p0<f
  11.             x0=x;p0=f;
  12.         end
  13.     end
  14. end
  15. x0,p0,toc
  16.  
  17. function [f,g]=mengte(x);
  18. f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)^2-8*x(1)-2*x(2)-3*x(3)-x(4)-2*x(5);
  19. g=[sum(x)-400
  20.     x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800
  21.     2*x(1)+x(2)+6*x(3)-200
  22.     x(3)+x(4)+5*x(5)-200];
  23. end
  24.  
  25. x0 =
  26.     46    98     1    99     3
  27. p0 =
  28.        50273

2.LINGO

  1. model:
  2. sets:
  3. row/1..4/:b;
  4. col/1..5/:c1,c2,x;
  5. link(row,col):a;
  6. endsets
  7. data:
  8. c1=1,1,3,4,2;
  9. c2=-8,-2,-3,-1,-2;
  10. a=1 1 1 1 1
  11. 1 2 2 1 6 
  12. 2 1 6 0 0 
  13. 0 0 1 1 5;
  14. b=400,800,200,200;
  15. enddata
  16. max=@sum(col:c1*x^2+c2*x);
  17. @for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
  18. @for(col:@gin(x));
  19. @for(col:@bnd(0,x,99));
  20. end

3.案例+灵敏度分析

（司P42）

  1. clc,clear,format long g,close all
  2. syms x1 x2
  3. f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2);
  4. f=simplify(f);
  5. f1=diff(f,x1);f2=diff(f,x2);
  6. [x10,x20]=solve(f1,f2);
  7. x10=round(double(x10));x20=round(double(x20));
  8. f0=subs(f,{x1,x2},{x10,x20});
  9. f0=double(f0);
  10. subplot(121),fmesh(f,[0,10000,0,10000]),title('')
  11. xlabel('$x_1$','Interpreter','latex')
  12. ylabel('$X_2$','Interpreter','latex')
  13. subplot(122),L=fcontour(f,[0,10000,0,10000]);
  14. contour(L.XData,L.YData,L.ZData,'ShowText','on')
  15. xlabel('$x_1$','Interpreter','latex')
  16. ylabel('$X_2$','Interpreter','latex','Rotation',0)
  17. p1=339-0.01*x10-0.003*x20
  18. p2=399-0.004*x10-0.01*x20
  19. c=400000+195*x10+225*x20        %总支出
  20. rate=f0/c                        %利润
  21. x10
  22. x20

灵敏度分析

  1. clc,clear,format long g,close all
  2. syms x1 x2 a
  3. f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(400000+195*x1+225*x2);
  4. f=simplify(f);
  5. f1=diff(f,x1);f2=diff(f,x2);
  6. [x10,x20]=solve(f1,f2);
  7. subplot(121),fplot(x10,[0.002,0.02]),title('')              %x1关于a的曲线
  8. xlabel('$x_1$','Interpreter','latex')
  9. ylabel('$X_2$','Interpreter','latex','Rotation',0)
  10. subplot(122),fplot(x20,[0.002,0.02]),title('')              %x2关于a的曲线
  11. xlabel('$x_1$','Interpreter','latex')
  12. ylabel('$X_2$','Interpreter','latex','Rotation',0)
  13. dx1=diff(x10,a);dx10=subs(dx1,a,0.01);dx10=double(dx10);
  14. sx1a=dx10*0.01/4735;
  15. dx2=diff(x20,a);dx20=subs(dx2,a,0.01);dx20=double(dx20);
  16. sx2a=dx20*0.01/7043;
  17. F=subs(f,{x1,x2},{x10,x20});
  18. F=simplify(F);
  19. figure,fplot(F,[0.002,0.02]),title('')
  20. xlabel('$x_1$','Interpreter','latex')
  21. ylabel('$X_2$','Interpreter','latex','Rotation',0)
  22. Sya=-4735^2*0.01/553641;
  23. f3=subs(f,{x1,x2,a},{4735,7043,0.011});
  24. f3=double(f3);
  25. f4=subs(F,a,0.011)                                          %近似最优利润
  26. f4=double(f4)                                               %最优利润
  27. delta=(f4-f3)/f4                                            %利润的相对误差

四.多目标规划

1.多目标规划

max f1(x)=2x(1)+3x(2)
min f2(x)=x(1)+2x(2)

0.5x(1)+0.25x(2)<=8
0.2x(1)+0.2x(2)<=4
x(1)+x(2)>=10
x(1)>=0
x(2)>=0

【1】线性加权

min 0.5(-f1)+0.5(f2)
结果：
x(1)=7，x(2)=13

【2】理想点

分别求解得f1=-53，f2=10
新的目标函数 min f=[-2x(1)-3x(2)+53]^2+[x(1)+2x(2)-10]^2
结果：
x(1)=13.36，x(2)=5.28

【3】序贯

分别求解得f1=-53
加入条件-2x(1)-3x(2)=-53求f2
结果：
x(1)=7，x(2)=13 
  1. clc,clear,prob=optimproblem;
  2. x=optimvar('x',2,'LowerBound',0);
  3. c1=[-2,-3];c2=[1,2];
  4. a=[0.5,0.25;0.2,0.2;1,5;-1,-1];
  5. b=[8;4;72;-10];
  6. prob.Constraints.con1=a*x<=b
  7. obj1=0.5*c1*x+0.5*c2*x
  8. prob1=prob;prob1.Objective=obj1;
  9. [sol1,fval1]=solve(prob1),sx=sol1.x
  10. f1=-c1*sx,f2=c2*sx
  11.  
  12. prob21=prob;prob21.Objective=c1*x;
  13. [sol21,fval21]=solve(prob21),sx21=sol21.x
  14. prob22=prob;prob22.Objective=c2*x;
  15. [sol22,fval22]=solve(prob22),sx22=sol22.x
  16. prob23=prob;
  17. prob23.Objective=(c1*x-fval21)^2+(c2*x-fval22)^2;
  18. [sol23,fval23]=solve(prob23),sx23=sol23.x
  19.  
  20. prob3=prob;prob3.Objective=c2*x;
  21. prob3.Constraints.con2=c1*x==fval21;
  22. [sol3,fval3]=solve(prob3),sx3=sol3.x

2.目标规划

正负偏差变量

软硬约束

“尽量”

优先权

序贯算法：

根据优先级排序，将目标规划分解成一系列单目标规划，依次求解

  1. clc,clear
  2. x=optimvar('x',2,'LowerBound',0);
  3. dp=optimvar('dp',4,'LowerBound',0);
  4. dm=optimvar('dm',4,'LowerBound',0);
  5. p=optimproblem('ObjectiveSense','min');
  6. p.Constraints.con1=2*sum(x)<=12;
  7. con2=[200*x(1)+300*x(2)+dm(1)-dp(1)==1500
  8.     2*x(1)-x(2)+dm(2)-dp(2)==0
  9.     4*x(1)+dm(3)-dp(3)==16
  10.     5*x(2)+dm(4)-dp(4)==15];
  11. p.Constraints.con2=con2;
  12. goal=100000*ones(3,1);
  13. mobj=[dm(1);dp(2)+dm(2);3*dp(3)+3*dm(3)+dp(4)];
  14. for i=1:3
  15.     p.Constraints.con3=[mobj<=goal];
  16.     p.Objective=mobj(i);
  17.     [sx,fval]=solve(p);
  18.     fprintf('第%d级目标计算结果如下:\n',i)
  19.     fval,xx=sx.x,sdm=sx.dm,sdp=sx.dp
  20.     goal(i)=fval;
  21. end

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/我家小花儿/article/detail/938582
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号