数学建模 算法与应用 （python实现）第二章 整数规划

第二章 整数规划

整数规划一般指整数线性规划

1. 概述

整数规划的分类

整数线性规划一般可分为两类：

1. 完全整数规划
2. 混合整数规划

整数规划的特点

1. 原线性规划有最优解
   • 原线性规划最优解全是整数，最优解一致
   • 整数规划无可行解
   • 有可行解，但最优解值会变差
2. 整数规划最优解不能按照师叔最优解简单的取证获得

求解方法分类

1. 分枝界定法——可求纯或混合整数线性规划。
2. 割平面法——可求纯或混合整数线性规划。
3. 隐枚举法——可求“0-1”整数规划。
   • 过滤隐枚举法
   • 分枝隐枚举法
4. 匈牙利法——解决指派问题（“0-1”规划特殊情形）。
5. 蒙特卡洛法——求解各种类型规划。

2. 0-1 型线性规划

0-1规划的变量$x_j$仅取值0或1。

相互排斥的约束条件

有两种运输方式可供选择，但只能选一种：
$$5x_1+4x_2\le 24或7x_1+3x_2\le 45\text{\hspace{1em}(2.1)}$$
为解决问题，引入0-1变量$y$
$$y=\{\begin{array}{l}1,当采用船运方法时\\ 0,当采用车运方法时\end{array}\text{\hspace{1em}(2.2)}$$
则上述条件改写为
$$\{\begin{array}{l}5x_1+4x_2\le 24+yM,\\ 7x_1+3x_2\le 45+yM,\\ y=0或1\end{array}\text{\hspace{1em}(2.3)}$$
其中M为充分大的数。

关于固定费用的问题(Fixed Cost Problem)

在某要求成本最小化的问题中，往往需要投入固定成本，例如，选定的生产方式投资高，由于产量大，因而分配到每件产品的成本就降低，反之初期成本投入低，因而每件产品的生产成本就偏高。

• $j=1,2,3$分别表示三种方式：

• $x_j$ 表示采用第 j 方式的产量

• $c_j$表示采用第$j$种方式时每件产品的成本变动

• $k_j$表示第$j$种方式时的固定成本

暂不考虑其他约束条件，采用各种方式的总成本为：
$$P_j=\{\begin{array}{l}{k}_j+c_j{x}_j,当x_j>0时\\ 0,当x_j=0时,j=1,2,3\end{array}\text{\hspace{1em}(2.4)}$$
为了统一讨论，引入0-1变量$y_j$，令
$$y_j=\{\begin{array}{l}1，当采用第i种生产方式时\\ 0，当不采用该生产方式时\end{array}\text{\hspace{1em}(2.5)}$$
于是目标函数转化为
$$\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)\text{\hspace{1em}(2.6)}$$
公式$(2.4)$可表示为下述三个线性约束条件
$$y_{j}m\le x_j\le y_{j}M\text{\hspace{1em}(2.7)}$$
式中$m$为充分小二点正常数；$M$为充分大的正常数

指派问题的数学模型

分配$n$人去做$n$项工作，每人做且仅作一项工作，若分配第$i$人去做第$j$项工作，需花费$C_{ij}$单位时间，如何分配时工人花费的总时间最少

引入0-1变量：
$$x_{ij}=\{\begin{array}{l}1,第i人做第j项工作\\ 0，第i人不做第j项工作\end{array}\text{\hspace{1em}(2.8)}$$
上述问的数学模型为
$$\begin{gathered} \min \sum _{i=1}^n\sum _{j=1}^n{c}_{ij}xij, \\ s.t.\{\begin{array}{l}{\sum }_{j=1}^n{x}_{ij}=1,i=1,...,n,\\ \\ {\sum }_{i=1}^n{x}_{ij}=1,j=1,...,n,\\ \\ x_{ij}=0或1,i,j=1,...,n.\end{array} \end{gathered}$$
求出最佳的x矩阵即可得到最优解，可使用匈牙利算法，拍卖算法等算法。

匈牙利算法python代码：

from scipy.optimize import linear_sum_assignment
  
cost =np.array([[4,1,3],[2,0,5],[3,2,2]])
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]
# [1 0 2] 
# [1 2 2] 
# 5

3. 蒙特卡洛法（随机取样法)

蒙特卡洛法也成为计算机随机模拟法，使用该方法必须使用计算机生成相关分布的随机数，

例：

$y=x^2,y=12-x与x轴在第一象限围成一个曲边三角形。设计一个随机试验，求该图形面积的近似值。$

import random

n = 10000000

x_min, x_max = 0.0, 12.0
y_min, y_max = 0.0, 9.0

total = 0.0

for i in range(n):
	x = random.uniform(x_min, x_max)
	y = random.uniform(y_min, y_max)
	if (y < x * x and x <= 3) or (y < 12 - x and x >= 3):
		total += 1

print(12*9*total/10**7)

# 输出：
49.511304

输出应在49.5附近，由于时随机模拟，因此每次输出结果都不完全相同

在一定计算量的情况下，蒙特卡洛法完全可以的出一个满意解。

例:

对于非线性整数规划，可使用蒙特卡洛法进行估算：
$$\begin{gathered} \max z=x_1^2+x_2^2+3x_3^2+4x_4^2+2x_5^2-8x_1-2x_2-3x_3-x_4-2x_5 \\ s.t.=\{\begin{array}{l}0\le x_i\le 99,\\ x_1+x_2+x_3+x_4+x_5\le 400,\\ x_1+2x_2+2x_3+x_4+6x_5\le 800,\\ 2x_1+x_2+6x_3\le 200,\\ x_3+x_4+5x_5\le 500.\end{array}\text{\hspace{1em}(2.9)} \end{gathered}$$

由于python下np随机耗时较大，使用C++

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define REP(i, lim) for(int i=0;i<lim;++i)
const int inf = 0x3f3f3f3f;

int x[10];

int check() {

    int sum = 0;
    REP(i, 5) sum+=x[i];

    if(sum>400) return 0;
    if(x[0]+2*x[1]+2*x[2]+x[3]+6*x[4] > 800) return 0;
    if(2*x[0]+x[1]+6*x[2]>200) return 0;
    if(x[2]+x[3]+5*x[4]>200) return 0;

    return 1;
}

void getRadom(int i) {
    srand((unsigned)time(NULL)*i);
    REP(i, 5)   x[i] = rand()%100;
    while(!check()) REP(i, 5) x[i] = rand()%100;
}

ll get_tot() {
    ll ans = x[0]*x[0] + x[1]*x[1] + 3*x[2]*x[2] + 4*x[3]*x[3] + 2*x[4]*x[4] -
            8*x[0]+2*x[1]-3*x[2]-x[3]-2*x[4];
    return ans;
}

int main()
{
    int lim = 2e7;
    ll tot = -inf;

    REP(i, lim) {
        getRadom(i);
        ll tmp = get_tot();
        tot = tot>tmp ? tot:tmp;
    }

    cout<<tot<<endl;

    return 0;
}
/*
输出：
51549
*/    

python代码:

import numpy as np


def check(x):
	if x.sum() > 400:
		return False
	if x[0]+2*x[1]+2*x[2]+x[3]+6*x[4] > 800:
		return False
	if 2*x[0]+x[1]+6*x[2]>200:
		return False
	if x[2]+x[3]+5*x[3]>200:
		return False

	return True


def get_radom():
	x = np.random.randint(100, size=5)
	while not check(x):
		x = get_radom()
	return x


lim = 10**6
ans = -1

for i in range(lim):
	num = get_radom()
	ans = max(ans, num.all())
	if i % 10000 == 0:
		print(i)

print('ans=' + ans)

太慢了没跑出结果(￣、￣)

4. 整数线性规划的计算机求解：

指派问题

上文提到了匈牙利算法的python实现及用法

混合整数规划

scipy包中没有相关的函数调用，因此只能使用分支界定法求解：

#FOR MY LITTLE CUTE BABY You
import math
from scipy.optimize import linprog
import sys
 
#数据的输入，数据的输入部分因为要使用linprog函数，所以输入的格式要注意，在下面我把linprog函数的作用以及函数的参数输入方式写下来。
 
c= [-3,-4,-1,-2]
A = [[2,3,1,3]]
b=[4]
x0_bounds=(0,1)
x1_bounds=(0,1)
x2_bounds=(0,1)
x3_bounds=(0,1)
 
#定义一个可以利用分支界定法的函数
def mylove(c,A,b,bounds=(x0_bounds,x1_bounds,x2_bounds,x3_bounds),t=1.0E-12):
    res = linprog(c, A_ub=A, b_ub=b,bounds=(x0_bounds,x1_bounds,x2_bounds,x3_bounds))
    #因为会有超出界限的取法，所以定义一个最大值，向下继续循环
    if (type(res.x) is float):
        bestX=-1*[sys.maxsize]*len(c)
    else:
        bestX=res.x
    bestVal=sum([x*y for x,y in zip(c, bestX)])
    if all(((x-math.floor(x))<t or (x-math.ceil(x))>t)  for x in bestX):
        return (bestVal,bestX)
    #如果不满足全部0-1之间整数输出的条件，就证明解里面有小数，这个时候需要把该小数化为0或者1
    else:
        ind = [i for i, x in enumerate(bestX) if (x-math.floor(x))>t and (x-math.ceil(x))<t][0]#如果当前数字是小数的判断
        newCon1=[0]*len(A[0])
        newCon2=[0]*len(A[0])
        newCon1[ind]=-1#把该变量添加条件为-xi<-math.ceil(bestX[ind]) 即Xi会取1
        newCon2[ind]=1#把该变量添加限制条件为xi<math.floor(bestX[ind]) 即xi会取0
        newA1=A.copy()
        newA2=A.copy()
        newA1.append(newCon1)
        newA2.append(newCon2)
        newB1=b.copy()
        newB2=b.copy()
        newB1.append(-math.ceil(bestX[ind]))
        newB2.append(math.floor(bestX[ind]))
        r1=mylove(c,newA1,newB1,bounds=(x0_bounds,x1_bounds,x2_bounds,x3_bounds))#两个分支不断往下面搜索，搜索全部为整数，最小的一个解答
        r2=mylove(c,newA2,newB2,bounds=(x0_bounds,x1_bounds,x2_bounds,x3_bounds))
        if r1[0]<r2[0]:
            return r1
        else:
            return r2
print(mylove(c,A,b,bounds=(x0_bounds,x1_bounds,x2_bounds,x3_bounds)))