这么详细的数学建模编程语言教程，确定不收藏吗？（二）

文案：shuohaodeBoltne

排版：随心390

hello，大家好。各位可点击左下方阅读原文，访问公众号官方店铺。谨防上当受骗，感谢各位支持！

​​这么详细的数学建模编程语言教程，确定不收藏吗？（一）​​这篇推文你们都看了吗？想必环境都配置好了吧？

啥？没配置好还在刷手手机机机（手动狗头）！！！进入主题

参考：知乎https://zhuanlan.zhihu.com/p/259849151

这只是进阶的一个问题，模型不多介绍、算法你们一看就懂

决策变量：

 车辆是否经过节点到达节点 

 车辆是否经过节点 

 车辆在访问节点后的剩余容量 

目标函数： 

约束条件：

01 | MATLAB+YALMIP+Gurobi

MATLAB基于Yalmip调用solver(Gurobi、Cplex)是非常容易上手的一门语言，对于运筹优化建模爱好者小菜一碟，直接上代码。

clear
clc
tic
V=5;  %车数量
N=13; %总节点数量
C=6; %单车容量
demands=[0,1.2,1.7,1.5,1.4,1.7,1.4,1.2,1.9,1.8,1.6,1.7,1.1]; %需求量
x=[81.5,87,75,85,89,77,76,87,73,77,73,91,92];
y=[41.5,37,53,52,41,58,45,53,38,38,31,47,44];
axis=[x' y']; %城市坐标
for i=1:N
    for j=1:N
        Dij(i,j)= sqrt((axis(i,1)-axis(j,1))^2+(axis(i,2)-axis(j,2))^2);
    end
end
%% 决策变量
xijk=binvar(N,N,V,'full');%i、j节点之间是否由第k辆车进行配送
yik=binvar(N,V,'full'); %第k辆车是否经过i节点
zik=sdpvar(N,V,'full'); %车辆k在访问i节点后，车子的剩余容量
%% 目标函数
obj=0;
for i=1:N
    for j=1:N
        for k=1:V
            obj=obj+Dij(i,j)*xijk(i,j,k);
        end
    end
end
%% 约束条件
Con=[ ];
for i = 2:N 
    Con = [Con; sum(yik(i,:)) == 1];
end


Con = [Con; sum(yik(1,:)) <= V];


for k = 1:V
    Con = [Con; sum(demands(:).*yik(:,k)) <= C];
end


for i = 1:N
    for j = 1:N
        for k = 1:V
            if i == j
                Con = [Con; xijk(i,j,k) == 0];
            end
        end
    end
end


for i = 1:N
    for k = 1:V
        Con = [Con; sum(xijk(i,:,k)) == sum(xijk(:,i,k))];
    end
end


for j = 1:N
    for k = 1:V
        Con = [Con; sum(xijk(:,j,k)) == yik(j,k)];
    end
end


for i = 1:N
    for k = 1:V
        Con = [Con; sum(xijk(i,:,k)) == yik(i,k)];
    end
end


for i = 2:N
    for j = 2:N
        for k = 1:V
            if i ~= j
                if demands(i)+demands(j) <= C
                    Con = [Con; zik(i,k)-zik(j,k)+C*xijk(i,j,k) <= C-demands(i)];
                end
            end
        end
    end
end


for i = 2:N
    for k = 1:V
        Con = [Con; zik(i,k) <= C];
        Con = [Con; zik(i,k) >= demands(i)];
    end
end
%% 求解
ops = sdpsettings( 'solver','gurobi');
sol = solvesdp(Con,obj,ops);
obj = double(obj);
xijk = double(xijk);
yik = double(yik);
zik = value(zik);
optimize(Con,obj)

注意：刚到图书馆时不要刷手机，要去逛逛YALMIP论坛。

02 | Python+Gurobi

机器学习，深度学习，大数据这么hot，确定不来了解一下罒ω罒。还是直接上代码，你们肯定都看过官方文档了，我就不在这班门弄斧了，小小小丢人！

from gurobipy import *
from scipy.spatial.distance import pdist, squareform
v = 5
n = 13
c = 6
demands = [0,1.2,1.7,1.5,1.4,1.7,1.4,1.2,1.9,1.8,1.6,1.7,1.1]
xy = [(81.5,41.5), (87,37), (75,53), (85,52), (89,41), (77,58), (76,45), (87,53), (73,38), (77,38), (73,31), (91,47), (92,44)]
distances = squareform(pdist(xy))
try:
    m = Model('cvrp')
    # 决策变量
    I = range(1, n+1)
    J = range(1, n+1)
    K = range(1, v+1)
    x = m.addVars(I, J, K, vtype = GRB.BINARY, name = "x")
    y = m.addVars(I, K, vtype = GRB.BINARY, name = 'y')
    z = m.addVars(I, K, vtype = GRB.CONTINUOUS, name = "z")
    # 更新变量
    m.update()
    # 目标函数
    m.setObjective(sum(distances[i-1,j-1]*x[i,j,k] for i in I for j in J for k in K), GRB.MINIMIZE)
    # 约束条件
    m.addConstrs(sum(y[i, k] for k in K) == 1 for i in I[1 :])
    m.addConstr(sum(y[1, k] for k in K) <= v)
    m.addConstrs(sum(demands[i-1]*y[i, k] for i in I) <= c for k in K)
    m.addConstrs(x[i,j, k] == 0 for i in I for j in J for k in K if i == j)
    m.addConstrs(sum(x[i,j,k] for j in J) == sum(x[j,i,k] for j in I) for i in I for k in K)
    m.addConstrs(sum(x[i,j,k] for i in I) == y[j,k] for j in J for k in K)
    m.addConstrs(sum(x[i,j,k] for j in J) == y[i,k] for i in I for k in K)
    m.addConstrs(z[i,k]-z[j,k]+c*x[i,j,k] <= c-demands[i-1] for i in I[1:] for j in J[1:] for k in K if i != j and demands[i-1]+demands[j-1] <= c)
    m.addConstrs(z[i,k] <= c for i in I[1:] for k in K)
    m.addConstrs(z[i,k] >= demands[i-1] for i in I[1:] for k in K)


    #求解
    m.optimize()
except GurobiError:
print('error reported')

注意注意：刚到图书馆时不要刷手机，取而代之的是去看官方文档。

03 | Julia+Gurobi

Julia是近年兴起的语言，在国内知名度很小，但官网称其有 C的速度，确定不来了解一下？PS：如果你在学习Julia之前学过任意一门语言，上手是不成问题的。废话不多说，见代码。

using JuMP
using Gurobi
v=5
n=13
c=6
demands=[0 1.2 1.7 1.5 1.4 1.7 1.4 1.2 1.9 1.8 1.6 1.7 1.1]
xy=[81.5 41.5;87 37;75 53;85 52;89 41;77 58;76 45;87 53;73 38;77 38;73 31;91 47;92 44]
distances = zeros(n,n)
for i in 1:n
    for j in 1:n
        global distances[i,j]=sqrt((xy[i,1]-xy[j,1])^2+(xy[i,2]-xy[j,2])^2)  
    end
end
distances
# M = Model(with_optimizer(CPLEX.Optimizer))
M = Model(Gurobi.Optimizer)
@variable(M,x[1:n,1:n,1:v],binary=true)
@variable(M,y[1:n,1:v],binary=true)
@variable(M,u[1:n,1:v]>=0)
for i in 2:n
   @constraint(M,sum(y[i,k] for k in 1:v) == 1)
end
@constraint(M,sum(y[1,k] for k in 1:v) <= v)
for k in 1:v
   @constraint(M,sum(demands[i]*y[i,k] for i in 1:n) <= c) 
end
for i in 1:n
    for k in 1:v
        @constraint(M,x[i,i, k] == 0)
        @constraint(M,sum(x[i,j,k] for j in 1:n) == sum(x[j,i,k] for j in 1:n))
    end
end
for j in 1:n
    for k in 1:v
        @constraint(M,sum(x[i,j,k] for i in 1:n) == y[j,k])
    end
end
for i in 1:n
    for k in 1:v
        @constraint(M,sum(x[i,j,k] for j in 1:n) == y[i,k])
    end
end
for i in 2:n
    for j in 2:n
        for k in 1:v
            if i != j
                if demands[i]+demands[j] <= c
                    @constraint(M,u[i,k]-u[j,k]+c*x[i,j,k] <= c-demands[i])
                end
            end
        end
    end
end
for i in 2:n
    for k in 1:v
        @constraint(M,u[i,k] <= c)
        @constraint(M,u[i,k] >= demands[i])
    end
end
obj = sum(distances[i,j]*x[i,j,k] for i in 1:n for j in 1:n for k in 1:v)
@objective(M,Min,obj)
# println(M)
optimize!(M)
objective_value.(M)
value.(x)
value.(y)
value.(u)

注意注意注意：刚到图书馆时不要刷手机，要去逛逛Julia官方论坛、官方文档。