数学建模(2)整数规划模型

整数规划模型^[1]

例 1

\[ \begin{gather} \max z=\sum_{i=1}^6x_i \\ \text{s.t.} \begin{cases} x_1+x_2\geq40 \\ x_2+x_3\geq50 \\ x_3+x_4\geq45 \\ x_4+x_5\geq55 \\ x_5+x_6\geq30 \\ x_6+x_1\geq35 \\ x_i\geq0,x_i\in\mathbb{R},\ i=1,2,\dots,6 \end{cases} \end{gather} \]

clear;
prob = optimproblem;
x = optimvar('x', 6, 'Type', 'integer', 'LowerBound', 0);
prob.Objective = sum(x);
con = optimconstr(6);

a = [40, 50, 45, 55, 30, 35];

for i = 1:5
    con(i) = x(i) + x(i + 1) >= a(i);
end

con(6) = x(6) + x(1) >= a(6);

prob.Constraints.con = con;
[sol, fval, flag, out] = solve(prob);
sol.x
fval

---

例 2

\[ \begin{gather} \max z=\sum_{i=1}^4\sum_{j=1}^5c_{ij}x_{ij} \\ \text{s.t.} \begin{cases} \sum_{i=1}^4x_{ij}=1, & j=1,2,\dots,5 \\ \sum_{j=1}^5x_{ij}\leq2, & i=1,2,3,4 \\ x_i=0\text{ or }1, & i=1,2,3,4,j=1,2,\dots,5 \end{cases} \end{gather} \]

clear;
c = load('data2.txt');
prob = optimproblem;
x = optimvar('x', 4, 5, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
prob.Objective = sum(sum(c .* x));
prob.Constraints.con1 = sum(x) == 1; %每列的和，返回行向量
prob.Constraints.con2 = sum(x, 2) <= 2; %每行的和，返回列向量
[sol, fval, flag, out] = solve(prob);
sol.x
fval

---

例 3：网点覆盖

\[ \begin{gather} \max z=\sum_{i=1}^{10}x_{i} \\ \text{s.t.} \begin{cases} \sum_{i=1}^{10}y_{ij}\geq1, & j=1,2,\dots,10 \\ d_{ij}y_{ij}\leq10x_i, & i,j=1,2,\dots,10 \\ \sum_{j=1}^{10}y_{ij}\leq5, & i=1,2,\dots,10 \\ x_i\geq y_{ij}, & i,j=1,2,\dots,10 \\ x_i=y_{ii}, & i=1,2,\dots,10 \\ \sum_{j=1}^5x_{ij}\leq2, & i=1,2,3,4 \\ x_i,y_{ij}=0\text{ or }1, & i,j=1,2,\dots,10 \end{cases} \end{gather} \]

clear;
a = load('data3.txt');
d = dist(a) %求两两列向量之间的距离
n = 10;
prob = optimproblem;
x = optimvar('x', n, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1)
y = optimvar('y', n, n, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
prob.Objective = sum(x);
con1 = [1 <= sum(y)'; sum(y, 2) <= 5];
con2 = optimconstr(2 * n ^ 2);
con3 = optimconstr(n);
k = 0;

for i = 1:n

    for j = 1:n
        k = k + 1;
        con2(2 * k - 1) = d(i, j) * y(i, j) <= 10 * x(i);
        con2(2 * k) = y(i, j) <= x(i);
    end

end

for i = 1:n
    con3(i) = x(i) == y(i, i);
end

prob.Constraints.con1 = con1;
prob.Constraints.con2 = con2;
prob.Constraints.con3 = con3;
[sol, fval, flag, out] = solve(prob);
xx = sol.x
yy = sol.y
fval

---

蒙特卡洛法

求 \(y=x^2\)，\(y=-x+12\) 与 \(x\) 轴在第一象限围成的面积。

clear;
n = 10 ^ 7;
x = unifrnd(0, 12, [1, n]);
y = unifrnd(0, 9, [1, n]);
freq = sum(y < x .^ 2 & x <= 3) + sum(y < 12 - x & x >= 3);
area_appr = 12 * 9 * freq / n

非线性整数规划：

clear;
%rng('shuffle')  %根据当前时间为随机数生成器提供种子
rng(0) %进行一致性比较，每次产生的随机数是一样的
p0 = 0;
n = 10 ^ 6;
tic %计时开始

for i = 1:n
    x = randi([0, 99], 1, 5); %产生一行五列的区间[0,99]上的随机整数
    [f, g] = MonteCarlo(x);

    if all(g <= 0)

        if p0 < f
            x0 = x;
            p0 = f; %记录下当前较好的解
        end

    end

end

x0, p0, toc %计时结束

function [f, g] = MonteCarlo(x); %定义目标函数和约束条件
    f = x(1) ^ 2 + x(2) ^ 2 + 3 * x(3) ^ 2 + 4 * x(4) ^ 2 + 2 * x(5) - 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

---

2005 年电工杯数学建模竞赛 B 题(整数规划，TSP 问题)

\[ \begin{gather} \max z=\sum_{i=1}^{15}\sum_{j=1}^{15}w_{ij}x_{ij} \\ \text{s.t.} \begin{cases} \sum_{j=1}^{15}x_{ij}=1, & i=1,2,\dots,15 \\ \sum_{i=1}^{15}x_{ij}=1, & j=1,2,\dots,15 \\ u_{i}-u_{j}+15x_{ij}\leq14, & i=1,2,\dots,15,j=2,3,\dots,15 \\ u_1=0, 1\leq u_\leq14, & i=2,3,\dots,15 \\ x_{ij}=0\text{ or }1, & i,j=1,2,\dots,15 \end{cases} \end{gather} \]

clear;
a = importdata('data6.xlsx');
a(isnan(a)) = 0; %把不确定值NaN替换为0
M = 10 ^ 7;
w = ones(14) * M;

for i = 1:14

    for j = 1:14
        if i ~= j, w(i, j) = sum(a(:, i) .* a(:, j)); end
    end

end

n = 15;
w(n, n) = M;
prob = optimproblem;
x = optimvar('x', n, n, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
u = optimvar('u', n, 'LowerBound', 0) %序号变量
prob.Objective = sum(sum(w .* x));
con1 = [sum(x, 2) == 1; sum(x)' == 1; u(1) == 0];
con2 = optimconstr(n * (n - 1));
k = 0;

for i = 1:n

    for j = 2:n
        k = k + 1;
        con2(k) = u(i) - u(j) + n * x(i, j) <= n - 1;
    end

end

prob.Constraints.con1 = con1;
prob.Constraints.con2 = con2;
prob.Constraints.con3 = [1 <= u(2:end)'; u(2:end)' <= 14];
[sol, fval, flag, out] = solve(prob);
[i, j] = find(sol.x);
ij = [i'; j']

---

例 4(Python)

为了便于对模型进行求解与分析，假设有 4 件事 4 个人去做，各变量对应的数据假设如表 1。

人	任务 A	任务 B	任务 C	任务 D
甲	25	29	31	42
乙	39	38	26	20
丙	34	27	28	40
丁	24	42	36	23

from scipy.optimize import linear_sum_assignment
import numpy as np

cost = np.array([[25, 29, 31, 42], [39, 38, 26, 20],
                 [34, 27, 28, 40], [24, 42, 36, 23]])
row_ind, col_ind = linear_sum_assignment(cost)
print(row_ind)  # 开销矩阵对应的行索引
print(col_ind)  # 对应行索引的最优指派的列索引
print(cost[row_ind, col_ind])  # 提取每个行索引的最优指派列索引所在的元素，形成数组
print(cost[row_ind, col_ind].sum())  # 数组求和

• 输出：

[0 1 2 3]
[0 2 1 3]
[25 26 27 23]
101

参考文献

1. 司守奎，孙玺菁. 数学建模算法与应用(第 3 版). ↩︎