cplex入门系列（四）--- 运输问题

一、运输问题理论介绍
在经济建设中，经常碰到大宗物资调运问题。如煤炭、钢材、木材、粮食等物资，在全国有若干生产基地，根据已有的交通网，应该如何指定调运方案，将这些物资运到各消费地点，而总运费要最小。这个问题可以用以下数学语言描述。
已知有 $m$ 个生产地点 $A_i,i=1,2,…,m$ 。可供应某种物资，其供应量（产量）分别为 $a_i,i=1,2,…,m$ ，有 $n$ 个销地 $B_j,j=1,2,…,n$ ，其需要量分别为 $b_j,j=1,2,…,n$ ，从 $A_i$ 到 $B_j$ 运输单位物资的运价（单价）为 $c_{ij}$ ，这些数据可以汇总于产销平衡表和单位运价表中。

产地	销地	产量
	1 2 … n
1		$a_1$
2		$a_2$
$⋮$		$⋮$
m		$a_m$
销量	$b_1 b_2…b_n$

产销平衡表

产地	销地
	1 2 … n
1	$c_{11}$ $c_{12}$ $…$ $c_{1n}$
2	$c_{21}$ $c_{22}$ $…$ $c_{2n}$
$⋮$	$...$
m	$c_{m1}$ $c_{m2}$ $…$ $c_{mn}$

单位运价表
若用 $x_ij$ 表示从 $A_i$ 到 $B_j$ 的运量，那么在产销平衡的条件下，要求得总运费最小得调运方案，可求解以下数学模型：
$min⁡ z=∑_{i=1}^m∑_{j=1}^nc_{ij}x_{ij}$
$\begin{equation} \left\{ \begin{array}{r1} \sum_{i=1}^mx_{ij}=b_j, j=1,2,\dots,n \\ \sum_{i=1}^nx_{ij}=a_i, i=1,2,\dots,m\\ x_{ij}{\geq}0 \end{array} \right. \end{equation}$
这就是运输问题得数学模型。它包含m×n个变量，(m+n)个约束方程。其系数矩阵得结构比较松散，且特殊。
该系数矩阵中对应于变量 $x_{ij}$ 的系数向量 $P_{ij}$ ，其分量中除第 $i$ 个和第 $m+j$ 个为1以外，其余的都为零。
$P_{ij}=(0…1…0…1…0)^T=e_i+e_{m+j}$
对产销平衡得运输问题，由于有以下关系存在
$∑_{j=1}^nb_j=∑_{j=1}^n(∑_{i=1}^mx_{ij})=∑_{i=1}^m(∑_{j=1}^nx_{ij})=∑_{i=1}^ma_i$
所以模型最多只有 $m+n-1$ 个独立约束方程。即稀疏矩阵得秩 $≤m+n−1$ 。由于有以上特征，所以求解运输问题时，可用比较简单得计算方法，习惯上称为表上作业法。
对产销平衡得运输问题，一定存在可行解。又因为所有变量有界，因而存在最优解。
二、cplex解决运输问题
归根到底运输问题还是一个线性规划问题，和cplex入门系列（二）中的LP问题并没有什么不同，但是由于运输问题是一个特殊的线性规划问题，即其一定存在可行解，且其结构特殊，故可以使用表上作业法进行求解。但是在cplex中，我们还是会使用求解一般线性规划的求解方式来求解问题。
2.1 案例一
假设有三个工厂Sunnyvale、Dublin、Bangkok供应四个地区Amarillo、Teaneck、Chicago和Sioux Falls的货物。各地区需求量及从各工厂到各地区送单位货物（柜）的运价表如图所示。试求出总的运费最节省的货物运输方案。

工厂/需求地区	Amarillo	Teaneck	Chicago	Sioux Falls	产量（柜）
Sunnyvale	250	420	380	280	45
Dublin	1280	990	1440	1520	120
Bangkok	1550	1420	1660	1730	95
需求量（柜）	80	78	47	55

数学建模：这是一个产销平衡的运输问题，对应数学模型为：
$min⁡ z=∑_{i=1}^3∑_{j=1}^4c_{ij}x_{ij}$
$\begin{equation} \left\{ \begin{array}{r1} \sum_{i=1}^3x_{ij}=b_j, j=1,2,3,4 \\ \sum_{i=1}^4x_{ij}=a_i, i=1,2,3\\ x_{ij}{\geq}0 \end{array} \right. \end{equation}$
用单位运价表和产销平衡表的数据带入上述数学模型
即有
$min⁡ z=250x_{11}+420x_{12}+380x_{13}+280x_{14}+1280x_{21}+990x_{22}+1440x_{23}+1520x_{24}+1550x_{31}+1420x_{32}+1660x_{33}+1730x_{34}$
$\begin{equation} \left\{ \begin{array}{r1} x_{11}+x_{12}+x_{13}+x_{14}=45 \\ x_{21}+x_{22}+x_{23}+x_{24}=120\\ x_{31}+x_{32}+x_{33}+x_{34}=95\\ x_{11}+x_{21}+x_{31}=80\\ x_{12}+x_{22}+x_{32}=78\\ x_{13}+x_{23}+x_{33}=47\\ x_{14}+x_{24}+x_{34}=55\\ \end{array} \right. \end{equation}$
该问题是一个纯粹的线性规划问题，即可以利润cplex的LP求解器来求解。

产/销地类定义（TransportationNode ）

package com.opreation.research;

/**
 * 这个用于定义运输节点
 */

public class TransportationNode {
    final String NodeName; // 节点名称
    final int quantity; // 产/销量
    final boolean isSource; // true 为产地 false 为销地
    final int id; // 下标

    public TransportationNode(String NodeName,int quantity,boolean isSource,int id) {
        this.NodeName = NodeName;
        this.quantity = quantity;
        this.isSource = isSource;
        this.id = id;
    }

    public boolean isSource() {
        return this.isSource;
    }
}

运价/距定义（TransportationRelation ）

package com.opreation.research;

/**
 * @Description: 定义运价、运距
 * @Author
 * @Date 2022/8/5 10:01
 */
public class TransportationRelation {
    final int [] distance ;
    public TransportationRelation(int[] distance) {
        this.distance = distance;
    }
}

cplex数据建模

package com.opreation.research;

import ilog.concert.IloException;
import ilog.concert.IloNumExpr;
import ilog.concert.IloNumVar;
import ilog.cplex.IloCplex;

import java.util.List;

/**
 * @Description: 定义运输模型
 * @Author 
 * @Date 2022/8/5 10:13
 */
public class TransportationModel {
    IloCplex cplex ;
    double objectiveValue;
    IloNumVar x[];
    int numOfSource = 0; // 产地数量
    int numOfSink = 0; // 销地数量
    /**
     * 假设有
     */
    public void bulidModel(List<TransportationNode> transportationNodeList, TransportationRelation transportationRelation)throws IloException {
        this.cplex = new IloCplex();
        for(TransportationNode tsNode : transportationNodeList){
            if (tsNode.isSource()) {
                numOfSource++;
            }else {
                numOfSink++;
            }
        }

        System.out.println(numOfSource);
        System.out.println(numOfSink);

        CreatDecisionVariable();
        CreatObjFunc(transportationRelation);
        CreatConstraints(transportationNodeList);
        Solve(transportationNodeList);
    }
    /**
     * 构建决策变量
     */
    public void CreatDecisionVariable() throws IloException{
        x = new IloNumVar[numOfSource*numOfSink];
        for (int i = 0; i < numOfSource; i++) {
            for (int j = 0; j < numOfSink; j++) {
                x[i*numOfSink+j] = cplex.numVar(0, Double.MAX_VALUE,"x"+(i+1)+(j+1));
                int a = i*numOfSink+j;
                System.out.println("x"+a+":" + "x"+(i+1)+(j+1));
            }
        }

    }
    /**
     * 构建目标函数，最小化运输距离（运价）
     */
    public void CreatObjFunc(TransportationRelation transportationRelation) throws IloException{
        cplex.addMinimize(cplex.scalProd(x,transportationRelation.distance));

    }
    /**
     * 构建约束条件
     */
    public void CreatConstraints(List<TransportationNode> transportationNodeList) throws IloException {
        for(TransportationNode tsNode : transportationNodeList){
            IloNumExpr left = cplex.numExpr();
            if (tsNode.isSource()) {
                // 产地节点，约束条件为，流出产地的供应量和产地决策变量相等
                for(int i = 0;i<numOfSink;i++) {
                    left = cplex.sum(left,cplex.prod(1,x[tsNode.id*numOfSink+i]));
                }
                cplex.addEq(left,tsNode.quantity);
            }else {
                // 销地节点，约束条件为，流入销地的供应量与销地决策变量相等
                for(int i = 0;i<numOfSource;i++) {
                    left = cplex.sum(left,cplex.prod(1,x[tsNode.id+i*numOfSink]));
                }
                cplex.addEq(left,tsNode.quantity);
            }
        }
    }
    /**
     * 线性规划模型求解问题
     */
    public void Solve(List<TransportationNode> transportationNodeList) throws IloException {
        cplex.setOut(null);
        if (cplex.solve()) {
            cplex.exportModel("transportationModel.lp");
            objectiveValue = cplex.getObjValue();
            System.out.println("最小运费为：" + objectiveValue);
            for (int i = 0; i < numOfSource; i++) {
                for (int j = 0; j < numOfSink; j++) {
                    System.out.print(transportationNodeList.get(i).NodeName+" to "+transportationNodeList.get(numOfSource+j).NodeName+" : ");
                    System.out.println(cplex.getValue(x[i*numOfSink+j]));
                }
            }
        }else{
            System.out.println("无解");
        }
    }
}

主函数（传入参数调用求解器）

package com.opreation.research;

import ilog.concert.IloException;

import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

/**
 * @Description: 传入参数，求解模型
 * @Author 
 * @Date 2022/8/5 10:56
 */
public class Main {
    public static void main(String[] args) throws IOException, IloException {
        // 产销数据节点
        TransportationNode transportationNode1 = new TransportationNode("Sunnyvale", 45, true, 0);
        TransportationNode transportationNode2 = new TransportationNode("Dublin", 120, true, 1);
        TransportationNode transportationNode3 = new TransportationNode("Bangkok", 95, true, 2);
        TransportationNode transportationNode4 = new TransportationNode("Amarillo", 80, false, 0);
        TransportationNode transportationNode5 = new TransportationNode("Teaneck", 78, false, 1);
        TransportationNode transportationNode6 = new TransportationNode("Chicago", 47, false, 2);
        TransportationNode transportationNode7 = new TransportationNode("Sioux Falls", 55, false, 3);
        List<TransportationNode> transportationNodeList = new ArrayList<TransportationNode>();
        transportationNodeList.add(transportationNode1);
        transportationNodeList.add(transportationNode2);
        transportationNodeList.add(transportationNode3);
        transportationNodeList.add(transportationNode4);
        transportationNodeList.add(transportationNode5);
        transportationNodeList.add(transportationNode6);
        transportationNodeList.add(transportationNode7);

        // 产销地运价/运距
        int [] dis = new int[]{250, 420, 380, 280, 1280, 990, 1440, 1520, 1550, 1420, 1660, 1730};
        TransportationRelation transportationRelation = new TransportationRelation(dis);

        TransportationModel model = new TransportationModel();
        model.bulidModel(transportationNodeList, transportationRelation);
    }
}

运行结果

最小运费为：297800.0
Sunnyvale to Amarillo : 0.0
Sunnyvale to Teaneck : 0.0
Sunnyvale to Chicago : 0.0
Sunnyvale to Sioux Falls : 45.0
Dublin to Amarillo : 42.0
Dublin to Teaneck : 78.0
Dublin to Chicago : 0.0
Dublin to Sioux Falls : 0.0
Bangkok to Amarillo : 38.0
Bangkok to Teaneck : 0.0
Bangkok to Chicago : 47.0
Bangkok to Sioux Falls : 10.0

打开transportationModel.lp可以看到模型结构，和我们构建的数学模型一致：

\ENCODING=ISO-8859-1
\Problem name: ilog.cplex

Minimize
 obj: 250 x11 + 420 x12 + 380 x13 + 280 x14 + 1280 x21 + 990 x22 + 1440 x23
      + 1520 x24 + 1550 x31 + 1420 x32 + 1660 x33 + 1730 x34
Subject To
 c1: x11 + x12 + x13 + x14  = 45
 c2: x21 + x22 + x23 + x24  = 120
 c3: x31 + x32 + x33 + x34  = 95
 c4: x11 + x21 + x31  = 80
 c5: x12 + x22 + x32  = 78
 c6: x13 + x23 + x33  = 47
 c7: x14 + x24 + x34  = 55
End

2.2 案例2
该案例为《运筹学》清华大学第四版P107 例4-1题目，具体应用场景请参考教材。
某公司经销售甲产品。它下设三个加工厂。每日的产量分别是： $A_1$ 为7吨， $A_2$ 为4吨， $A_3$ 为9吨。该公司把这些产品分别运往四个销售点。各销售点每日销量为： $B_1$ 为3吨， $B_2$ 为6吨， $B_3$ 为5吨， $B_4$ 为6吨。已知从各工厂到各销售点的单位产品的运价如下表所示。问该公司应该如何调运产品，在满足各销点的需求量的前提下，使总运费为最少。

加工厂/销地	$B_1$	$B_2$	$B_3$	$B_4$
$A_1$	3	11	3	10
$A_2$	1	9	2	8
$A_3$	7	4	10	5

加工厂/销地	$B_1$	$B_2$	$B_3$	$B_4$	产量
$A_1$					7
$A_2$					4
$A_3$					9
销量	3	6	5	6

数学建模：这是一个产销平衡的运输问题，对应数学模型为：
$min⁡ z=∑_{i=1}^3∑_{j=1}^4c_{ij}x_{ij}$
$\begin{equation} \left\{ \begin{array}{r1} \sum_{i=1}^3x_{ij}=b_j, j=1,2,3,4 \\ \sum_{i=1}^4x_{ij}=a_i, i=1,2,3\\ x_{ij}{\geq}0 \end{array} \right. \end{equation}$
用单位运价表和产销平衡表的数据带入上述数学模型
即有
$min⁡ z=3x_{11}+11x_{12}+3x_{13}+10x_{14}+1x_{21}+9x_{22}+2x_{23}+8x_{24}+7x_{31}+4x_{32}+10x_{33}+5x_{34}$
$\begin{equation} \left\{ \begin{array}{r1} x_{11}+x_{12}+x_{13}+x_{14}=7\\ x_{21}+x_{22}+x_{23}+x_{24}=4\\ x_{31}+x_{32}+x_{33}+x_{34}=9\\ x_{11}+x_{21}+x_{31}=3\\ x_{12}+x_{22}+x_{32}=6\\ x_{13}+x_{23}+x_{33}=5\\ x_{14}+x_{24}+x_{34}=6\\ \end{array} \right. \end{equation}$
该问题是一个纯粹的线性规划问题，即可以利用cplex的LP求解器来求解。
产/销地类定义（TransportationNode ）

package com.opreation.research;

/**
 * 这个用于定义运输节点
 */

public class TransportationNode {
    final String NodeName; // 节点名称
    final int quantity; // 产/销量
    final boolean isSource; // true 为产地 false 为销地
    final int id; // 下标

    public TransportationNode(String NodeName,int quantity,boolean isSource,int id) {
        this.NodeName = NodeName;
        this.quantity = quantity;
        this.isSource = isSource;
        this.id = id;
    }

    public boolean isSource() {
        return this.isSource;
    }
}