数学建模——线性规划类-Toy模板网

这篇具有很好参考价值的文章主要介绍了数学建模——线性规划类。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

一.线性规划

1.算法

【1】通用代码

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

例如：

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

matlab

（1）基于求解器

%    LP1_1_1.m

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    %最大值


x =
    6.4286
    0.5714
         0
y =
   14.5714
% 说明：在x1 x2 x3 = 6.4286 0.5714 0 的情况下，maxz = 14.5714

（2）基于问题

con中根据符号分类

%    LP1_1_2.m

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;
x=sol.x,y=fval

python

#    线性规划模型.py

from scipy.optimize import linprog
import numpy as np

c=np.array([-2,-3,5])
a=[[-2,5,-1],[1,3,1]]
b=[[-10],[12]]
aeq=[[1,1,1]]
beq=[7]
lb=[0,0,0]
ub=[None]*len(c)
bound=tuple(zip(lb,ub))
res=linprog(c,a,b,aeq,beq,bound,method='simplex',options={"disp":True})
print("最优解：",res.x)
print("目标函数最小值：",res.fun)
print("目标函数最大值：",-res.fun)
print(res)

最优解： [6.42857143 0.57142857 0.        ]
目标函数最小值： -14.571428571428571
目标函数最大值： 14.571428571428571
     con: array([0.])
     fun: -14.571428571428571
 message: 'Optimization terminated successfully.'
     nit: 3
   slack: array([0.        , 3.85714286])
  status: 0
 success: True
       x: array([6.42857143, 0.57142857, 0.        ])

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

（1）绝对值

%    LP_abs.m

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


x =
    -2
     0
     0
     0
minz =
    2.0000

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

（见风投案例模型二）

2.风投案例

【0】题目描述

【1】模型一

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

（1）matlab

%    LP_fengtou_1.m

clc,clear
t=0;hold on
while t<0.05
    c=[-0.05,-0.27,-0.19,-0.185,-0.185];
    A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];
    b=t*ones(4,1);
    Aeq=[1,1.01,1.02,1.045,1.065];beq=1;
    lb=zeros(5,1);
    [x,Q]=linprog(c,A,b,Aeq,beq,lb);
    Q=-Q;
    plot(t,Q,'*k');
    t=t+0.001;
end
xlabel('t'),ylabel('Q')

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

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

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

（2）python

#    风投案例.py

from scipy.optimize import linprog
import numpy as np
from numpy import ones,diag,c_,zeros
import matplotlib.pyplot as plt

plt.rc('font',size=16)
c=np.array([-0.05,-0.27,-0.19,-0.185,-0.185])
A=c_[zeros(4),diag([0.025,0.015,0.055,0.026])]
Aeq=[[1,1.01,1.02,1.045,1.065]]
beq=[1]
a=0
aa=[]
ss=[]
while a<0.05:
    b=ones(4)*a
    res = linprog(c, A, b, Aeq, beq)
    x=res.x
    Q=-res.fun
    aa.append(a)
    ss.append(Q)
#输出转折点
    if a==0.006:
        print("x=", x)
        print("Q=", Q)
    a=a+0.001
plt.plot(aa,ss,'r*')
plt.xlabel('$a$')
plt.ylabel('$Q$',rotation=90)
plt.savefig('模型一.png',dpi=500)
plt.show()

x= [0.         0.24       0.4        0.10909091 0.22122066]
Q= 0.20190763977806234（万元）

【2】模型二

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

import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp

plt.rc('font',size=16)
x=cp.Variable(6,pos=True)
obj=cp.Minimize(x[5])
a1=np.array([0.025,0.015,0.055,0.026])
a2=np.array([0.05,0.27,0.19,0.185,0.185])
a3=np.array([1,1.01,1.02,1.045,1.065])
k=0.05
kk=[]
ss=[]
while k<0.27:
    con=[cp.multiply(a1,x[1:5])-x[5]<=0,
         0.05*x[0]+0.27*x[1]+0.19*x[2]+0.185*x[3]+0.185*x[4]>=k,
         x[0]+1.01*x[1]+1.02*x[2]+1.045*x[3]+1.065*x[4]==1]
    prob=cp.Problem(obj,con)
    prob.solve(solver='CVXOPT')
#输出转折点
    if abs(k-0.21)<0.0005:
        print("x=",x.value)
        print("风险=",prob.value) 
    kk.append(k)
    ss.append(prob.value)
    k=k+0.005
plt.plot(kk,ss,'r*')
plt.xlabel('$k$')
plt.ylabel('$R$',rotation=90)
plt.savefig('模型二.png',dpi=500)
plt.show()

x= [2.09147087e-08 3.08874349e-01 5.14790372e-01 1.40396702e-01
 1.52452147e-02 7.72185738e-03]
风险= 0.00772185738406626

【3】模型三

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

（1）matlab

clc,clear
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('收益')
    X
end
V,Q,format

w=（权重）
    [0.766,    0.767,    0.810,    0.811,    0.824,    0.825,    0.962,    0.963,    1.0];
V =（风险）
  247.5248   92.2509   92.2509   78.4929   78.4929   59.3960   59.3960         0         0
Q =（收益）
   1.0e+03 *
    2.6733    2.1648    2.1648    2.1060    2.1060    2.0162    2.0162    0.5000    0.5000
X =
   1.0e+04 *

         0    0.9901         0         0         0    0.0248
         0    0.3690    0.6150         0         0    0.0092
         0    0.3690    0.6150         0         0    0.0092
         0    0.3140    0.5233    0.1427         0    0.0078
         0    0.3140    0.5233    0.1427         0    0.0078
         0    0.2376    0.3960    0.1080    0.2284    0.0059（转折点处的X值）
         0    0.2376    0.3960    0.1080    0.2284    0.0059
    1.0000         0         0         0         0         0
    1.0000         0         0         0         0         0

（2）python

from scipy.optimize import linprog
import matplotlib.pyplot as plt

plt.rc('font',size=16)
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]]
b=[0,0,0,0]
Aeq=[[1,1.01,1.02,1.045,1.065,0]]
beq=[1]
bound=((0,None),(0,None),(0,None),(0,None),(0,None),(0,None))
s=0;ss=[];aa=[]
while s<=1:
    c=[-(1-s)*0.05,-(1-s)*0.27,-(1-s)*0.19,-(1-s)*0.185,-(1-s)*0.185,s]
    res=linprog(c,A,b,Aeq,beq)
    ss.append(s);aa.append(-res.fun)
    s=s+0.01
plt.plot(ss,aa,'r.')
plt.show()

import matplotlib.pyplot as plt
import cvxpy as cp
import numpy as np

plt.rc('font',size=16)
x=cp.Variable(6)
a1=np.array([0.025,0.015,0.055,0.026])
a2=np.array([0.05,0.27,0.19,0.185,0.185])
a3=np.array([1,1.01,1.02,1.045,1.065])
con=[cp.multiply(a1,x[1:5])-x[5]<=0,a3@x[:-1]==1,x>=0]
s=0;ss=[];aa=[]
while s<=1:
    obj=cp.Minimize(s*x[-1]-(1-s)*a2@x[:-1])
    prob=cp.Problem(obj,con)
    prob.solve(solver='GLPK_MI')
    ss.append(s);aa.append(-prob.value)
    s=s+0.01
plt.plot(ss,aa,'r.')
plt.show()

二.整数规划

0.分类

1.通用代码

%    int_LP_1.m

clc,clear,prob=optimproblem;
x=optimvar('x',6,'Type','integer','LowerBound',0);
prob.Objective=sum(x);
% 创建一个由空优化约束组成的 6×1 数组，使用 constr 初始化用于创建约束表达式的循环
coon=optimconstr(6);
a=[35,40,50,45,55,30];
con(1)=x(1)+x(6)>=35;
for i=1:5
    con(i+1)=x(i)+x(i+1)>=a(i+1);
end
prob.Constraints.con=con;
[sol,fval,flag]=solve(prob)
x=sol.x

x =
    35
     5
    45
     0
    55
     0

intlinprog函数

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

%    int_LP_intlinprog.m

c = [-40;-90];
A = [9,7;7,20];
b = [56;70];
[x,fval] = intlinprog(c,[1 2],A,b,[],[],zeros(2,1));
x,maxz=-fval

x =
    4.0000
    2.0000
maxz =
   340

2.分支定界法（剪枝）

尝试一维单调

%    int_branchbunch.m
clear global;
clear;
clc;

global result; % 存储所有整数解
global lowerBound; % 下界
global upperBound; % 上界
global count; % 用以判断是否为第一次分支

count = 1;

f = [-40, -90];
A = [8, 7;7, 20;];  
b = [56; 70];
Aeq = [];
beq = [];
lbnd = [0; 0];
ubnd = [inf; inf];

% f = [-3 2 -5];
% A = [1 2 -1;1 4 1;1 1 0;0 4 1];
% b = [2;4;3;6];
% Aeq = [];
% beq = [];
% lbnd=[0 0 0];
% ubnd=[1 1 1];

BinTree = createBinTreeNode({f, A, b, Aeq, beq, lbnd, ubnd});
if ~isempty(result)
    [fval,flag]=min(result(:,end)); % result中每一行对应一个整数解及对应的函数值
    Result=result(flag,:);
    disp('该整数规划问题的最优解为：');
    disp(Result(1,1:end-1));
    disp('该整数规划问题的最优值为：');
    disp(Result(1,end));
else
    disp('该整数规划问题无可行解');
end

function BinTree = createBinTreeNode(binTreeNode)

global result;
global lowerBound;
global upperBound;
global count;

if isempty(binTreeNode)
    return;
else
    BinTree{1} = binTreeNode;
    BinTree{2} = [];
    BinTree{3} = [];
    
    [x, fval, exitflag] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
        binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, binTreeNode{7});
    if count == 1
%         upperBound = 0; % 初始下界为空
        lowerBound = fval;
        count = 2;
    end
    
    if ~IsInRange(fval)
        return;
    end

    if exitflag == 1
        flag = [];
        % 寻找非整数解分量
        for i = 1 : length(x)
            if round(x(i)) ~= x(i)
                flag = i;
                break;
            end
        end
        % 分支
        if ~isempty(flag)
            lowerBound = max([lowerBound; fval]);
            OutputLowerAndUpperBounds();
            lbnd = binTreeNode{6};
            ubnd = binTreeNode{7};
            lbnd(flag, 1) = ceil(x(flag, 1)); % 朝正无穷四舍五入
            ubnd(flag, 1) = floor(x(flag, 1));
            
            % 进行比较，优先选择目标函数较大的进行分支
            [~, fval1] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
        binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd);
            [~, fval2] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
        binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7});
            if fval1 < fval2                
                % 创建左子树          
                BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});

                % 创建右子树
                BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
            else
                % 创建右子树
                BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}});
                % 创建左子树          
                BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ...
            binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd});
            end
        else
            upperBound = min([upperBound; fval]);
            OutputLowerAndUpperBounds();
            result = [result; [x', fval]];
            return;
        end
    else
        % 剪枝
        return;
    end  
end
end

function y = IsInRange(fval)
    global lowerBound;
    global upperBound;

    if isempty(upperBound) & fval >= lowerBound
        y = 1;
    else if  (fval >= lowerBound & fval <= upperBound)
        y = 1;
    else
        y = 0;
    end
end
end

function y = OutputLowerAndUpperBounds()

global lowerBound;
global upperBound;

disp("此时下界、上界分别为");
disp(lowerBound);
if isempty(upperBound)
    disp('  无上界')
else
    disp(upperBound);
end

end

该整数规划问题的最优解为：
     4     2
该整数规划问题的最优值为：
  -340

3.割平面法

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

function  [intx,intf] = DividePlane(A,c,b,baseVector)
%功能：用割平面法求解整数规划
%调用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
%其中, A：约束矩阵；
%      c：目标函数系数向量；
%      b：约束右端向量；
%      baseVector：初始基向量；
%      intx：目标函数取最值时的自变量值；
%      intf：目标函数的最值；
sz = size(A);
nVia = sz(2);%获取有多少决策变量
n = sz(1);%获取有多少约束条件
xx = 1:nVia;

if length(baseVector) ~= n
    disp('基变量的个数要与约束矩阵的行数相等！');
    mx = NaN;
    mf = NaN;
    return;
end
 
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
 
%首先用单纯形法求出最优解
while 1   
    [maxs,ind] = max(sigma);
%--------------------用单纯形法求最优解--------------------------------------
    if maxs <= 0   %当检验数均小于0时，求得最优解。      
        vr = find(c~=0 ,1,'last');
        for l=1:vr
            ele = find(baseVector == l,1);
            if(isempty(ele))
                mx(l) = 0;
            else
                mx(l)=xb(ele);
            end
        end
        if max(abs(round(mx) - mx))<1.0e-7  %判断最优解是否为整数解，如果是整数解。
            intx = mx;
            intf = mx*c;
            return;
        else  %如果最优解不是整数解时，构建切割方程
            sz = size(A);
            sr = sz(1);
            sc = sz(2);
            [max_x, index_x] = max(abs(round(mx) - mx));
            [isB, num] = find(index_x == baseVector);
            fi = xb(num) - floor(xb(num));
            for i=1:(index_x-1)
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            for i=(index_x+1):sc
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            
            Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
            A = [A zeros(sr,1);-Atmp(1,:) 1];
            xb = [xb;-fi];
            baseVector = [baseVector sc+1];
            sigma = [sigma 0];
         
            %-------------------对偶单纯形法的迭代过程----------------------
            while 1
                %----------------------------------------------------------
                if xb >= 0    %判断如果右端向量均大于0，求得最优解
                    if max(abs(round(xb) - xb))<1.0e-7   
%如果用对偶单纯形法求得了整数解，则返回最优整数解
                        vr = find(c~=0 ,1,'last');
                        for l=1:vr
                            ele = find(baseVector == l,1);
                            if(isempty(ele))
                                mx_1(l) = 0;
                            else
                                mx_1(l)=xb(ele);
                            end
                        end
                        intx = mx_1;
                        intf = mx_1*c;
                        return;
                    else   %如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程
                        sz = size(A);
                        sr = sz(1);
                        sc = sz(2);
                        [max_x, index_x] = max(abs(round(mx_1) - mx_1));
                        [isB, num] = find(index_x == baseVector);
                        fi = xb(num) - floor(xb(num));
                        for i=1:(index_x-1)
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        for i=(index_x+1):sc
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        Atmp(1,index_x) = 0;  %下一次对偶单纯形迭代的初始表格
                        A = [A zeros(sr,1);-Atmp(1,:) 1];
                        xb = [xb;-fi];
                        baseVector = [baseVector sc+1];
                        sigma = [sigma 0];
                        continue;
                    end
                else   %如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
                    minb_1 = inf;
                    chagB_1 = inf;
                    sA = size(A);
                    [br,idb] = min(xb);
                    for j=1:sA(2)
                        if A(idb,j)<0
                            bm = sigma(j)/A(idb,j);
                            if bm<minb_1
                                minb_1 = bm;
                                chagB_1 = j;
                            end
                        end
                    end
                    sigma = sigma -A(idb,:)*minb_1;
                    xb(idb) = xb(idb)/A(idb,chagB_1);
                    A(idb,:) = A(idb,:)/A(idb,chagB_1);
                    for i =1:sA(1)
                        if i ~= idb
                            xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
                            A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
                        end
                    end
                    baseVector(idb) = chagB_1;
                end
              %------------------------------------------------------------
            end 
            %--------------------对偶单纯形法的迭代过程---------------------    
        end     
    else     %如果检验数有不小于0的，则进行单纯形算法的迭代过程
        minb = inf;
        chagB = inf;
        for j=1:n
            if A(j,ind)>0
                bz = xb(j)/A(j,ind);
                if bz<minb
                    minb = bz;
                    chagB = j;
                end
            end
        end
        sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
        xb(chagB) = xb(chagB)/A(chagB,ind);
        A(chagB,:) = A(chagB,:)/A(chagB,ind);
        for i =1:n
            if i ~= chagB
                xb(i) = xb(i)-A(i,ind)*xb(chagB);
                A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
            end
        end
        baseVector(chagB) = ind;
    end
    M = M + 1;
    if (M == 1000000)
        disp('找不到最优解！');
        mx = NaN; 
        minf = NaN;
        return;
    end
end

%    int_LP_divideplane.m

c = [-1;-1]; % 不要加松弛变量
A = [-1 1 1 0;3 1 0 1]; % 加上松弛变量
b = [1;4];
[x,fval] = DividePlane(A,c,b,[3 4]); % 松弛变量 3 4
x,maxz=fval

x =
    1.0000    1.0000
maxz =
    -2

4.匈牙利算法

【1】

【2】

%    int_XiongYaLi.m

c = [6 7 11 2;4 5 9 8;3 1 10 4;5 9 8 2];
a = zeros(8,16);
for i = 1:4
a(i,(i-1)*4+1:4*i) = 1;
a(4+i,i:4:16) = 1;
end
b = ones(8,1);
[x,y] = linprog(c,[],[],a,b,zeros(16,1),ones(16,1));
X = reshape(x,4,4)
opt = y

% X =
% 
%      0     0     0     1
%      1     0     0     0
%      0     1     0     0
%      0     0     1     0
% 
% 
% opt =
% 
%     15

A=[6 7 11 2;
    4 5 9 8;
    3 1 10 4;
    5 9 8 2];
% 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];
B=Hungary(A);
[~,b]=linear_assignment(A,B)

% b =
%        4              1              8              2 

function res=Hungary(N)
%输入的矩阵应N*N的
[a,~]=size(N);
%第一步每一行减去当前行最小值
for ii = 1:a
    N(ii,:)= N(ii,:)-min( N(ii,:));
end
%第二步每一列减去当前列最小值
for ii = 1:a
    N(:,ii)=  N(:,ii)-min( N(:,ii));
end
num=0;
while num~=a
    [num,N_min,del_hang,del_lie]=line_count(N);
    if num ~=a
        for ii=1:a
            if del_hang(ii)~=ii
                N(ii,:) =  N(ii,:)-N_min;
            end
            if del_lie(ii)==ii
            N(:,ii) =  N(:,ii)+N_min;
            end
        end
    else
        res=N;
    end
end

function [num,M_min,del_hang,del_lie]=line_count(M)
[a,~]=size(M);
num=0;
h=0;
del_hang=zeros(a,1);
del_lie=zeros(a,1);
for ii=1:a
    del=ii-h;
    [~,b]=size(find(M(del,:)==0));
    if   b>= 2
        M(del,:)=[];
        h=h+1;
        del_hang(ii)=ii;    %得到被覆盖的行数
        num=num+1;
    end
end
l=0;
for ii=1:a
    del=ii-l;
    [b,~]=size(find(M(:,del)==0));
    if  b >=1
        M(:,del)=[];
        l=l+1;
        del_lie(ii)=ii;    %得到被覆盖的列数
        num=num+1;
    end
end
M_min=min(min(M));
end
end

function [place,res]=linear_assignment(M,N)
%N是n维矩阵,N是经过Hungary处理的
%M是未处理前的
[a,~]=size(N);
x=0;
place=zeros(1,a);
res=zeros(1,a);
judge=zeros(1,a);
while find(N==0)
    for ii=1:a
    judge(ii)=length(find(N(ii,:)==0));
    end
    judge(find(judge==0))=[];
    if min(judge)==1
     for ii=1:a
        if length(find(N(ii,:)==0))==1     %先选出行中只有1个0
            x=x+1;
            place(x)=ii+(find(N(ii,:)==0)-1)*a; %得到矩阵中的位置
            h=find(N(ii,:)==0);
            N(ii,:)=1./zeros(1,a);
            N(:,h)=1./zeros(a,1);
        end
     end
    end
    
    for ii=1:a
    judge(ii)=length(find(N(ii,:)==0));
    end
    judge(find(judge==0))=[];
    
    if min(judge)==2
       x=x+1;
    q=find(N==0);
    place(x)=q(1);
    N(mod(q(1),a),:)=1./zeros(1,a);
    N(:,fix(q(1)/a)+1)=1./zeros(a,1);  
    end
end
[place,~]=sort(place);
for ii=1:length(place)
    res(ii)=M(place(ii));
end
end

5.0-1规划

【1】intlinprog函数

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

%    intlinprog_01.m

c = [-3 2 -5]';
A = [1 2 -1;1 4 1;1 1 0;0 4 1];
b = [2;4;3;6];
[x,fval] = intlinprog(c,[1 2 3],A,b,[],[],zeros(3,1),ones(3,1));
x,maxz = -fval

x =
     1
     0
     1
maxz =
     8

【2】指派问题

（1）纯0-1整数规划

c=np.array([[4,8,7,15,12],
            [7,9,17,14,10],
            [6,9,12,8,7],
            [6,7,14,6,10],
            [6,9,12,10,6]])
x=cp.Variable((5,5),integer=True)
obj=cp.Minimize(cp.sum(cp.multiply(c,x)))
con=[0<=x,x<=1,
     cp.sum(x,axis=0,keepdims=True)==1,
     cp.sum(x,axis=1,keepdims=True)==1]
prob=cp.Problem(obj,con)
prob.solve(solver='GLPK_MI')
print("最优值为：",prob.value)
print("最优解为：\n",x.value)

最优值为： 34.0
最优解为：
 [[0. 0. 1. 0. 0.]
 [0. 1. 0. 0. 0.]
 [1. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]

（2）广义指派

5*x(1)+4*x(2)<=24+y(1)M;        %M是充分大的数
7*x(1)+3*x(2)<=45+y(2)M;
y(1)+y(2)=1;

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

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

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));
y(i)*m<=x(i)<=y(i)*M;   %m充分小，M充分大
y(i)=0 or 1;

【3】旅行商问题(TSP)

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

三.非线性规划

算法

【1】二次规划

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

H为实对称矩阵

[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)

%    DLP_1.m

H=[4,-4;-4,8];
f=[-6;3];
a=[1,1;4,1];
b=[3;9];
[x,fval]=quadprog(H,f,a,b,[],[],zeros(2,1))

x =
    2.0854
    0.6585
fval =
   -5.5976

%    DLP_2.m

x0 = [1;1];
A=[1,1;4,1];
b=[3;9];
Aeq = [];
beq = [];
lb = [0;0];
[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb)
function f = fun(x)
    f = 2*x(1)-4*x(1)*x(2)+4*x(2)^2-6*x(1)-3*x(2)
end

x =
    2.0000
    1.0000
fval =
  -15.0000

【2】无约束非线性规划

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

%    NLP_nocon.m

clc,clear
f = @(x) x(1)^3-x(2)^3+3*x(1)^2+3*x(2)^2-9*x(1);
g = @(x) -f(x);
[m1,n1] = fminunc(f,[0,0])%求极小值
[m2,n2] = fminsearch(g,[0,0]);%求极大值
m2,n2=-n2

m1 =
    1.0000   -0.0000
n1 =
    -5
m2 =
   -3.0000    2.0000
n2 =
   31.0000

【3】有约束非线性规划

%    NLP_con_1.m
%[x,fval]=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub,'fun2')

[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')

function f = fun1(x);
f=sum(x.^2)+8;
end

function [g,h] = fun2(x)
g=[-x(1)^2+x(2)-x(3)^2
    x(1)+x(2)^2+x(3)^3-20];
h=[-x(1)-x(2)^2+2
    x(2)+2*x(3)^2-3];
end

x =
    0.5522
    1.2033
    0.9478
y =
   10.6511

%    NLP_con_2.m

clc,clear
x = optimvar('x',3,'LowerBound',0);
prob = optimproblem('Objective',sum(x.^2)+8);
con1 = [-x(1)^2+x(2)-x(3)^2 <= 0
    x(1)+x(2)^2+x(3)^3 <= 20];
con2 = [-x(1)-x(2)^2+2 == 0
    x(2)+2*x(3)^2 == 3];
prob.Constraints.con1 = con1;
prob.Constraints.con2 = con2;
x0.x = rand(3,1);
[s,f,flag,out] = solve(prob,x0);
s.x,f

（1）只有等式约束

拉格朗日法

（2）一般形式

罚函数法

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

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

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

%    NLP_TU.m

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

（2）K-T条件

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

若为凸规划，则充要条件

【5】蒙特卡洛法

%    MTKL.m

clc,clear
rng(0)          %rng('shuffle')
p0=0;n=10^6;tic
for i=1:n
    x=randi([0,99],1,5);
    [f,g]=mengte(x);
    if all(g<=0)
        if p0<f
            x0=x;p0=f;
        end
    end
end
x0,p0,toc

function [f,g]=mengte(x);
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);
g=[sum(x)-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];
end

x0 =
    46    98     1    99     3
p0 =
       50273

2.LINGO

model:
sets:
row/1..4/:b;
col/1..5/:c1,c2,x;
link(row,col):a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=1 1 1 1 1
1 2 2 1 6 
2 1 6 0 0 
0 0 1 1 5;
b=400,800,200,200;
enddata
max=@sum(col:c1*x^2+c2*x);
@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
@for(col:@gin(x));
@for(col:@bnd(0,x,99));
end

3.案例+灵敏度分析

（司P42）

clc,clear,format long g,close all
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);
[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),L=fcontour(f,[0,10000,0,10000]);
contour(L.XData,L.YData,L.ZData,'ShowText','on')
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
p1=339-0.01*x10-0.003*x20
p2=399-0.004*x10-0.01*x20
c=400000+195*x10+225*x20        %总支出
rate=f0/c                        %利润
x10
x20

灵敏度分析

clc,clear,format long g,close all
syms x1 x2 a
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);
[x10,x20]=solve(f1,f2);
subplot(121),fplot(x10,[0.002,0.02]),title('')              %x1关于a的曲线
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
subplot(122),fplot(x20,[0.002,0.02]),title('')              %x2关于a的曲线
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
dx1=diff(x10,a);dx10=subs(dx1,a,0.01);dx10=double(dx10);
sx1a=dx10*0.01/4735;
dx2=diff(x20,a);dx20=subs(dx2,a,0.01);dx20=double(dx20);
sx2a=dx20*0.01/7043;
F=subs(f,{x1,x2},{x10,x20});
F=simplify(F);
figure,fplot(F,[0.002,0.02]),title('')
xlabel('$x_1$','Interpreter','latex')
ylabel('$X_2$','Interpreter','latex','Rotation',0)
Sya=-4735^2*0.01/553641;
f3=subs(f,{x1,x2,a},{4735,7043,0.011});
f3=double(f3);
f4=subs(F,a,0.011)                                          %近似最优利润
f4=double(f4)                                               %最优利润
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

clc,clear,prob=optimproblem;
x=optimvar('x',2,'LowerBound',0);
c1=[-2,-3];c2=[1,2];
a=[0.5,0.25;0.2,0.2;1,5;-1,-1];
b=[8;4;72;-10];
prob.Constraints.con1=a*x<=b
obj1=0.5*c1*x+0.5*c2*x
prob1=prob;prob1.Objective=obj1;
[sol1,fval1]=solve(prob1),sx=sol1.x
f1=-c1*sx,f2=c2*sx

prob21=prob;prob21.Objective=c1*x;
[sol21,fval21]=solve(prob21),sx21=sol21.x
prob22=prob;prob22.Objective=c2*x;
[sol22,fval22]=solve(prob22),sx22=sol22.x
prob23=prob;
prob23.Objective=(c1*x-fval21)^2+(c2*x-fval22)^2;
[sol23,fval23]=solve(prob23),sx23=sol23.x

prob3=prob;prob3.Objective=c2*x;
prob3.Constraints.con2=c1*x==fval21;
[sol3,fval3]=solve(prob3),sx3=sol3.x

2.目标规划

正负偏差变量

软硬约束

“尽量”

优先权

序贯算法：

根据优先级排序，将目标规划分解成一系列单目标规划，依次求解文章来源地址https://www.toymoban.com/news/detail-649112.html

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