一、问题描述

选址问题是指在规划区域里选择一个或多个设施的位置，使得目标最优。

按照规划区域的结构划分，可以将选址问题分为：
1.连续选址问题：设施可以在给定范围的任意位置选址，设施的候选位置为无穷多。
2.离散选址问题：设施的候选位置是有限且较少的，实际中最常遇到这类问题。
3.网格选址问题：规划区域被划分为许多的小单元，每个设施占据其中有限个单元。

本博客将根据下面的例题，介绍连续选址问题的求解思路
（其中，问题1为运输问题，建议先阅读我的另一篇博客再来看本博客：【运筹优化】运输问题建模 + Java调用Cplex求解）

在这里插入图片描述

二、思路分析

问题（2）属于连续选址问题，它的数学模型和问题（1）的数学模型 $Q_1$ 类似，唯一不同之处在于，问题（2）中， $d_{ij}$ 也是一个决策变量。因此，它的数学模型是非线性的，无法直接通过Cplex直接求解。对于凸的无约束非线性模型，可以使用梯度下降法或牛顿法求出其最优解；对于其他非线性模型，通常使用启发式算法或群体智能算法寻求其较优解。由于不确定问题（2）模型是否为凸，因此本人采用了一个贪心启发式算法和一个群体智能算法(蜘蛛猴优化算法)对问题（2）进行求解。

三、解决方案

3.1 贪心启发式算法

首先介绍贪心启发式算法，为了更好地描述，给出下面的线性规划模型：

在这里插入图片描述
该模型的含义为：给定一个已知位置的料场和若干个工地的信息，确保将料场的储量全部运输到各个工地，并使得总运输成本最小。

然后介绍贪心启发式算法的基本思想：每次为一个料场 $i$ 定位置，基于每个工地的坐标和当前剩余需求量计算一个加权中心点作为当前料场 $i$ 的坐标，在确定了料场 $i$ 的位置之后，就可以将料场 $i$ 和当前工地的信息代入到线性的 $Q^i_2$ 模型中调用Cplex进行求解，得到料场 $i$ 的最小成本的运输方案后，再更新每个工地的剩余需求量，然后以同样的步骤为下一个料场确定位置并更新工地信息，直到每个工地的位置都被确定。确定了每个工地的位置之后，再求解 $Q_1$ 模型，即可得到当前料场坐标下最优的运输方案。

贪心启发式的Java代码实现如下：

public class Answer2 {
    /**
     * 料场对象
     */
    @AllArgsConstructor
    static class Stockyard {
        /**
         * x,y坐标和储量c
         */
        double x, y, c;
    }

    /**
     * 工地对象
     */
    @AllArgsConstructor
    static class ConstructionSite {
        /**
         * x,y坐标和需求d
         */
        double x, y, d;
    }

    private static double calcDistance(double x1, double y1, double x2, double y2) {
        return Math.sqrt(Math.pow(x1 - x2, 2) + Math.pow(y1 - y2, 2));
    }

    private static List<ConstructionSite> getInitConstructionSiteList() {
        List<ConstructionSite> constructionSiteList = new ArrayList<>();
        constructionSiteList.add(new ConstructionSite(1.25, 1.25, 3));
        constructionSiteList.add(new ConstructionSite(8.75, 0.75, 5));
        constructionSiteList.add(new ConstructionSite(0.5, 4.75, 4));
        constructionSiteList.add(new ConstructionSite(5.75, 5, 7));
        constructionSiteList.add(new ConstructionSite(3, 6.5, 6));
        constructionSiteList.add(new ConstructionSite(7.25, 7.75, 11));
        return constructionSiteList;
    }

    public static void main(String[] args) throws Exception {
        // 浮点型精度误差
        double EPS = 1e-06;
        // 料场列表
        List<Stockyard> stockyardList = new ArrayList<>();
        stockyardList.add(new Stockyard(0, 0, 20));
        stockyardList.add(new Stockyard(0, 0, 20));
        stockyardList.sort(new Comparator<Stockyard>() {
            @Override
            public int compare(Stockyard o1, Stockyard o2) {
                return -Double.compare(o1.c, o2.c);
            }
        });
        // 工地列表
        List<ConstructionSite> constructionSiteList = getInitConstructionSiteList();
        // 初始化距离矩阵
        double[][] distanceMatrix = new double[stockyardList.size()][constructionSiteList.size()];
        // 开始循环
        for (int i = 0; i < stockyardList.size(); i++) {
            // 计算总需求量
            double totalDemand = constructionSiteList.stream().map(constructionSite -> {
                return constructionSite.d;
            }).reduce(0d, Double::sum);
            // 计算加权中心点作为料场i的坐标
            stockyardList.get(i).x = 0d;
            stockyardList.get(i).y = 0d;
            for (ConstructionSite constructionSite : constructionSiteList) {
                stockyardList.get(i).x += (constructionSite.x * constructionSite.d / totalDemand);
                stockyardList.get(i).y += (constructionSite.y * constructionSite.d / totalDemand);
            }
            // 计算距离矩阵
            for (int j = 0; j < distanceMatrix[i].length; j++) {
                distanceMatrix[i][j] = calcDistance(stockyardList.get(i).x, stockyardList.get(i).y, constructionSiteList.get(j).x, constructionSiteList.get(j).y);
            }
            if (totalDemand <= stockyardList.get(i).c) {
                break;
            }
            // 然后强制料场i必须送够c_i的量，并计算最小成本的运输方案
            // 开始建模
            IloCplex cplex = new IloCplex();
            // 声明变量
            IloNumVar[] x = new IloNumVar[constructionSiteList.size()];
            for (int j = 0; j < x.length; j++) {
                x[j] = cplex.numVar(0, constructionSiteList.get(j).d);
            }
            // 构造约束1：不能超出每个工地的需求
            for (int j = 0; j < constructionSiteList.size(); j++) {
                cplex.addLe(x[j], constructionSiteList.get(j).d);
            }
            // 构造约束2：料场中的储量必须全部送出
            cplex.addEq(cplex.sum(x), stockyardList.get(i).c);
            // 声明目标函数
            cplex.addMinimize(cplex.scalProd(x, distanceMatrix[i]));
            // 配置cplex
            cplex.setOut(null);
            cplex.setWarning(null);
            cplex.setParam(IloCplex.DoubleParam.EpOpt, EPS);
            // 开始求解
            if (cplex.solve()) {
                // 更新每个工地的需求量
                for (int j = 0; j < x.length; j++) {
                    constructionSiteList.get(j).d -= cplex.getValue(x[j]);
                }
            } else {
                System.err.println("此题无解");
            }
            // 结束模型
            cplex.end();
        }
        for (int i = 0; i < stockyardList.size(); i++) {
            System.out.println("料场" + (i + 1) + "的坐标为: ( " + stockyardList.get(i).x + " , " + stockyardList.get(i).y + " )");
        }
        // 重新获取初始化的工地列表
        constructionSiteList = getInitConstructionSiteList();
        // 开始建模
        IloCplex cplex = new IloCplex();
        // 声明变量
        IloNumVar[][] x = new IloNumVar[stockyardList.size()][constructionSiteList.size()];
        for (int i = 0; i < x.length; i++) {
            for (int j = 0; j < x[i].length; j++) {
                x[i][j] = cplex.numVar(0, Math.min(stockyardList.get(i).c, constructionSiteList.get(j).d));
            }
        }
        // 构造约束1：必须满足每个工地的需求
        for (int j = 0; j < constructionSiteList.size(); j++) {
            IloLinearNumExpr expr = cplex.linearNumExpr();
            for (int i = 0; i < x.length; i++) {
                expr.addTerm(1, x[i][j]);
            }
            cplex.addEq(expr, constructionSiteList.get(j).d);
        }
        // 构造约束2：不能超过每个料场的储量
        for (int i = 0; i < stockyardList.size(); i++) {
            cplex.addLe(cplex.sum(x[i]), stockyardList.get(i).c);
        }
        // 声明目标函数
        IloLinearNumExpr target = cplex.linearNumExpr();
        for (int i = 0; i < x.length; i++) {
            for (int j = 0; j < x[i].length; j++) {
                target.addTerm(distanceMatrix[i][j], x[i][j]);
            }
        }
        cplex.addMinimize(target);
        // 配置cplex
        cplex.setOut(null);
        cplex.setWarning(null);
        cplex.setParam(IloCplex.DoubleParam.EpOpt, EPS);
        // 开始求解
        long s = System.currentTimeMillis();
        if (cplex.solve()) {
            System.out.println("最小吨千米数为: " + cplex.getObjValue());
            for (int i = 0; i < stockyardList.size(); i++) {
                for (int j = 0; j < x[i].length; j++) {
                    double xValue = cplex.getValue(x[i][j]);
                    if (xValue > EPS) {
                        System.out.println("料场" + (i + 1) + "向工地" + (j + 1) + "运输" + xValue + "吨水泥");
                    }
                }
            }
            System.out.println("求解用时: " + (System.currentTimeMillis() - s) / 1000d + " s");
        } else {
            System.err.println("此题无解");
        }
        // 结束模型
        cplex.end();
    }
}

运行结果如下：

料场1的坐标为: ( 5.208333333333333 , 5.159722222222222 )
料场2的坐标为: ( 4.90625 , 3.59375 )
最小吨千米数为: 114.45936860248577
料场1向工地4运输7.0吨水泥
料场1向工地5运输2.0吨水泥
料场1向工地6运输11.0吨水泥
料场2向工地1运输3.0吨水泥
料场2向工地2运输5.0吨水泥
料场2向工地3运输4.0吨水泥
料场2向工地5运输4.0吨水泥
求解用时: 0.0 s

3.2 群体智能算法（蜘蛛猴优化算法）

然后是使用了蜘蛛猴优化算法。每个蜘蛛猴的变量数组表示了每个料场的坐标。设置蜘蛛猴数量为50，迭代次数为60，局部领导者决策阶段的更新几率为0.1，局部领导者阶段的更新几率为0.8，最大组数为5，料场的坐标上下界分别为10和-10。

Java代码实现如下：

public class SMO_Solver {

    /**
     * 料场对象
     */
    @AllArgsConstructor
    static class Stockyard {
        /**
         * x,y坐标和储量c
         */
        double x, y, c;
    }

    /**
     * 工地对象
     */
    @AllArgsConstructor
    static class ConstructionSite {
        /**
         * x,y坐标和需求d
         */
        double x, y, d;
    }

    private double calcDistance(double x1, double y1, double x2, double y2) {
        return Math.sqrt(Math.pow(x1 - x2, 2) + Math.pow(y1 - y2, 2));
    }

    private static List<ConstructionSite> getInitConstructionSiteList() {
        List<ConstructionSite> constructionSiteList = new ArrayList<>();
        constructionSiteList.add(new ConstructionSite(1.25, 1.25, 3));
        constructionSiteList.add(new ConstructionSite(8.75, 0.75, 5));
        constructionSiteList.add(new ConstructionSite(0.5, 4.75, 4));
        constructionSiteList.add(new ConstructionSite(5.75, 5, 7));
        constructionSiteList.add(new ConstructionSite(3, 6.5, 6));
        constructionSiteList.add(new ConstructionSite(7.25, 7.75, 11));
        return constructionSiteList;
    }

    List<ConstructionSite> constructionSiteList;
    List<Stockyard> stockyardList;
    double EPS = 1e-06;

    public SMO_Solver(List<Stockyard> stockyardList, List<ConstructionSite> constructionSiteList) throws IloException {
        this.constructionSiteList = constructionSiteList;
        this.stockyardList = stockyardList;
    }

    // 蜘蛛猴对象
    class SpiderMonkey {
        // 当前蜘蛛猴坐标（自变量数组）
        double[] curVars;
        // 当前自变量对应的目标函数值
        double curObjValue;
        // 适应度（解决最小化问题，所以适应度为目标函数值的倒数）
        double fit;

        // 全参数构造
        public SpiderMonkey(double[] curVars, double curObjValue, double fit) {
            this.curVars = curVars;
            this.curObjValue = curObjValue;
            this.fit = fit;
        }
    }

    // 算法参数
    // 蜘蛛猴群
    List<SpiderMonkey[]> spiderMonkeyList = new ArrayList<>();
    // 局部领导者
    List<SpiderMonkey> localLeaderList = new ArrayList<>();
    // 最好的蜘蛛猴(全局领导者)
    SpiderMonkey bestSpiderMonkey;
    // 随机数对象
    Random random = new Random();
    // 最大迭代次数
    int maxGen = 60;
    // 蜘蛛猴数量
    int spiderMonkeyNum = 50;
    // 局部搜索次数(一般等于蜘蛛猴数量)
    int localSearchCount = spiderMonkeyNum;
    // 局部领导者决策阶段的更新几率
    double LLDP_PR = 0.1;
    // 局部领导者阶段的更新几率
    double LLP_PR = 0.8;
    // 变量维度数
    int varNum = 4;
    // 最大组数(一般要至少保证每组里有10个蜘蛛猴)
    int maxgroupNum = spiderMonkeyNum / 10;
    // 变量的上下界
    double[] ub = new double[]{10, 10, 10, 10};
    double[] lb = new double[]{-10, -10, -10, -10};
    // 局部计数器
    int[] localLimitCount = new int[]{0};
    // 停止条件
    int limitCnt = 30;
    // 全局计数器
    int globalLimitCount;
    // 记录迭代过程
    public double[][][] positionArr;
    // 记录迭代器的行数
    int curC = 0;
    // 是否开启贪心机制（只接受比当前解好的解）
    boolean greedy = true;

    // 求解主函数
    public void solve() {
        // 初始化蜘蛛猴种群
        initSpiderMonkeys();
        // 开始迭代
        for (int t = 0; t < maxGen; t++) {
            System.out.println(t + " : " + bestSpiderMonkey.curObjValue);
            // 局部领导者阶段（LLP：所有的蜘蛛猴都有机会更新自己）
            LLP();
            // 全局领导者阶段（GLP：轮盘赌，随机选取，偏向于对fit值大的蜘蛛猴进行更新）
            GLP();
            // 全局领导者学习阶段（如果全局领导者有更新，则globalLimitCount=0，否则globalLimitCount++）
            GLLP();
            // 局部领导者学习阶段（如果局部领导者有更新，则localLimitCount=0，否则localLimitCount++）
            LLLP();
            // 局部领导者决策阶段
            LLDP();
            // 全局领导者决策阶段
            GLDP();
        }
        // 输出最好的结果
        System.out.println("变量取值为：" + Arrays.toString(bestSpiderMonkey.curVars));
        System.out.println("最优解为：" + bestSpiderMonkey.curObjValue);
    }

    // 全局领导者决策阶段
    private void GLDP() {
        if (globalLimitCount >= limitCnt) {
            globalLimitCount = 0;
            if (spiderMonkeyList.size() < maxgroupNum) {
                // 分裂
                List<SpiderMonkey> tempList = new ArrayList<>();
                for (SpiderMonkey[] spiderMonkeys : spiderMonkeyList) {
                    tempList.addAll(Arrays.asList(spiderMonkeys));
                }
                tempList.sort(new Comparator<SpiderMonkey>() {
                    @Override
                    public int compare(SpiderMonkey o1, SpiderMonkey o2) {
                        return Double.compare(o2.fit, o1.fit);
                    }
                });
                //
                int groupNum = spiderMonkeyList.size() + 1;
                spiderMonkeyList = new ArrayList<>();
                int avgNum = spiderMonkeyNum / groupNum;
                for (int i = 0; i < groupNum - 1; i++) {
                    SpiderMonkey[] spiderMonkeys = new SpiderMonkey[avgNum];
                    for (int j = 0; j < avgNum; j++) {
                        spiderMonkeys[j] = copySpiderMonkey(tempList.remove(0));
                    }
                    spiderMonkeyList.add(spiderMonkeys);
                }
                spiderMonkeyList.add(tempList.toArray(new SpiderMonkey[0]));
                localLimitCount = new int[groupNum];
            } else {
                // 融合
                SpiderMonkey[] spiderMonkeys = new SpiderMonkey[spiderMonkeyNum];
                int i = 0;
                for (SpiderMonkey[] monkeys : spiderMonkeyList) {
                    for (SpiderMonkey monkey : monkeys) {
                        spiderMonkeys[i++] = copySpiderMonkey(monkey);
                    }
                }
                spiderMonkeyList = new ArrayList<>();
                spiderMonkeyList.add(spiderMonkeys);
                localLimitCount = new int[]{0};
            }
            // 更新局部领导者
            localLeaderList = new ArrayList<>();
            for (SpiderMonkey[] spiderMonkeys : spiderMonkeyList) {
                localLeaderList.add(copySpiderMonkey(spiderMonkeys[0]));
                int index = localLeaderList.size() - 1;
                for (int i = 1; i < spiderMonkeys.length; i++) {
                    if (localLeaderList.get(index).fit < spiderMonkeys[i].fit) {
                        localLeaderList.set(index, copySpiderMonkey(spiderMonkeys[i]));
                    }
                }
            }
        }
    }

    // 局部领导者决策阶段
    private void LLDP() {
        int c = 0;
        for (int i = 0; i < spiderMonkeyList.size(); i++) {
            SpiderMonkey[] spiderMonkeys = spiderMonkeyList.get(i);
            if (localLimitCount[i] < limitCnt) {
                for (int j = 0; j < spiderMonkeys.length; j++) {
                    SpiderMonkey tempSpiderMonkey = copySpiderMonkey(spiderMonkeys[j]);
                    for (int m = 0; m < varNum; m++) {
                        if (random.nextDouble() <= LLDP_PR) {
                            tempSpiderMonkey.curVars[m] = lb[m] + random.nextDouble() * (ub[m] - lb[m]);
                        } else {
                            double moveDist = random.nextDouble() * (bestSpiderMonkey.curVars[m] - tempSpiderMonkey.curVars[m]) + random.nextDouble() * (spiderMonkeys[random.nextInt(spiderMonkeys.length)].curVars[m] - tempSpiderMonkey.curVars[m]);
                            moveSpiderMonkey(tempSpiderMonkey, m, moveDist);
                        }
                    }
                    tempSpiderMonkey.curObjValue = getObjValue(tempSpiderMonkey.curVars);
                    tempSpiderMonkey.fit = 1 / tempSpiderMonkey.curObjValue;
                    if (greedy) {
                        if (tempSpiderMonkey.fit > spiderMonkeys[j].fit) {
                            spiderMonkeys[j] = tempSpiderMonkey;
                        }
                    } else {
                        spiderMonkeys[j] = tempSpiderMonkey;
                    }
                }
            }
            for (int j = 0; j < spiderMonkeys.length; j++) {
                for (int m = 0; m < spiderMonkeys[j].curVars.length; m++) {
                    positionArr[curC][c][m] = spiderMonkeys[j].curVars[m];
                }
                c++;
            }
        }
        curC++;
    }

    // 局部领导者学习阶段（如果局部领导者有更新，则localLimitCount=0，否则localLimitCount++）
    private void LLLP() {
        for (int i = 0; i < spiderMonkeyList.size(); i++) {
            boolean isUpdate = false;
            for (SpiderMonkey spiderMonkey : spiderMonkeyList.get(i)) {
                if (localLeaderList.get(i).fit < spiderMonkey.fit) {
                    localLeaderList.set(i, copySpiderMonkey(spiderMonkey));
                    isUpdate = true;
                }
            }
            if (isUpdate) {
                localLimitCount[i] = 0;
            } else {
                localLimitCount[i]++;
            }
        }
    }

    // 全局领导者学习阶段（如果全局领导者有更新，则globalLimitCount=0，否则globalLimitCount++）
    private void GLLP() {
        boolean isUpdate = false;
        for (int i = 0; i < spiderMonkeyList.size(); i++) {
            for (SpiderMonkey spiderMonkey : spiderMonkeyList.get(i)) {
                if (spiderMonkey.fit > bestSpiderMonkey.fit) {
                    bestSpiderMonkey = copySpiderMonkey(spiderMonkey);
                    isUpdate = true;
                }
            }
        }
        if (isUpdate) {
            globalLimitCount = 0;
        } else {
            globalLimitCount++;
        }
    }

    // 全局领导者阶段（GLP：轮盘赌，随机选取，偏向于对fit值大的蜘蛛猴进行更新）
    private void GLP() {
        int c = 0;
        for (int i = 0; i < spiderMonkeyList.size(); i++) {
            SpiderMonkey[] spiderMonkeys = spiderMonkeyList.get(i);
            // 计算fit总和
            double totalFit = 0;
            for (SpiderMonkey spiderMonkey : spiderMonkeys) {
                totalFit += spiderMonkey.fit;
            }
            // 轮盘赌的累计概率数组
            double[] p = new double[spiderMonkeys.length];
            for (int j = 0; j < p.length; j++) {
                p[j] = (spiderMonkeys[j].fit / totalFit) + (j == 0 ? 0 : p[j - 1]);
            }
            // 局部搜索
            for (int j = 0; j < localSearchCount; j++) {
                double r = random.nextDouble();
                for (int k = 0; k < p.length; k++) {
                    if (r <= p[k]) {
                        for (int m = 0; m < varNum; m++) {
                            double moveDist = random.nextDouble() * (bestSpiderMonkey.curVars[m] - spiderMonkeys[k].curVars[m]) + (random.nextDouble() - 0.5) * 2 * (spiderMonkeys[random.nextInt(spiderMonkeys.length)].curVars[m] - spiderMonkeys[k].curVars[m]);
                            moveSpiderMonkey(spiderMonkeys[k], m, moveDist);
                        }
                        spiderMonkeys[k].curObjValue = getObjValue(spiderMonkeys[k].curVars);
                        spiderMonkeys[k].fit = 1 / spiderMonkeys[k].curObjValue;
                        break;
                    }
                }
            }
            for (int j = 0; j < spiderMonkeys.length; j++) {
                for (int m = 0; m < spiderMonkeys[j].curVars.length; m++) {
                    positionArr[curC][c][m] = spiderMonkeys[j].curVars[m];
                }
                c++;
            }
            spiderMonkeyList.set(i, spiderMonkeys);
        }
        curC++;
    }

    // 局部领导者阶段（LLP：所有的蜘蛛猴都有机会更新自己）
    private void LLP() {
        int c = 0;
        for (int i = 0; i < spiderMonkeyList.size(); i++) {
            SpiderMonkey[] spiderMonkeys = spiderMonkeyList.get(i);
            SpiderMonkey localLeader = localLeaderList.get(i);
            for (int j = 0; j < spiderMonkeys.length; j++) {
                // 以一定几率更新自己
                if (random.nextDouble() <= LLP_PR) {
                    SpiderMonkey tempSpiderMonkey = copySpiderMonkey(spiderMonkeys[j]);
                    for (int m = 0; m < varNum; m++) {
                        double moveDist = random.nextDouble() * (localLeader.curVars[m] - tempSpiderMonkey.curVars[m]) + (random.nextDouble() - 0.5) * 2 * (spiderMonkeys[random.nextInt(spiderMonkeys.length)].curVars[m] - tempSpiderMonkey.curVars[m]);
                        moveSpiderMonkey(tempSpiderMonkey, m, moveDist);
                    }
                    tempSpiderMonkey.curObjValue = getObjValue(tempSpiderMonkey.curVars);
                    tempSpiderMonkey.fit = 1 / tempSpiderMonkey.curObjValue;
                    if (greedy) {
                        if (tempSpiderMonkey.fit > spiderMonkeys[j].fit) {
                            spiderMonkeys[j] = tempSpiderMonkey;
                        }
                    } else {
                        spiderMonkeys[j] = tempSpiderMonkey;
                    }
                }
                for (int m = 0; m < spiderMonkeys[j].curVars.length; m++) {
                    positionArr[curC][c][m] = spiderMonkeys[j].curVars[m];
                }
                c++;
            }
            spiderMonkeyList.set(i, spiderMonkeys);
        }
        curC++;
    }

    // 初始化蜘蛛猴种群
    private void initSpiderMonkeys() {
        positionArr = new double[3 * maxGen][spiderMonkeyNum][varNum];
        SpiderMonkey[] spiderMonkeys = new SpiderMonkey[spiderMonkeyNum];
        SpiderMonkey localLeader = null;
        for (int i = 0; i < spiderMonkeyNum; i++) {
            spiderMonkeys[i] = getRandomSpiderMonkey();
            if (i == 0) {
                bestSpiderMonkey = copySpiderMonkey(spiderMonkeys[0]);
                localLeader = copySpiderMonkey(spiderMonkeys[0]);
            } else {
                if (bestSpiderMonkey.fit < spiderMonkeys[i].fit) {
                    bestSpiderMonkey = copySpiderMonkey(spiderMonkeys[i]);
                    localLeader = copySpiderMonkey(spiderMonkeys[0]);
                }
            }
        }
        spiderMonkeyList.add(spiderMonkeys);
        localLeaderList.add(localLeader);
    }

    // 获取一个随机生成的蜘蛛猴
    SpiderMonkey getRandomSpiderMonkey() {
        double[] vars = new double[varNum];
        for (int j = 0; j < vars.length; j++) {
            vars[j] = lb[j] + random.nextDouble() * (ub[j] - lb[j]);
        }
        double objValue = getObjValue(vars);
        return new SpiderMonkey(vars.clone(), objValue, 1 / objValue);
    }

    // 控制spiderMonkey在第m个维度上移动n个距离
    public void moveSpiderMonkey(SpiderMonkey spiderMonkey, int m, double n) {
        // 移动
        spiderMonkey.curVars[m] += n;
        // 超出定义域的判断
        if (spiderMonkey.curVars[m] < lb[m]) {
            spiderMonkey.curVars[m] = lb[m];
        }
        if (spiderMonkey.curVars[m] > ub[m]) {
            spiderMonkey.curVars[m] = ub[m];
        }
    }

    /**
     * @param vars 自变量数组
     * @return 返回目标函数值
     */
    IloCplex cplex = new IloCplex();

    public double getObjValue(double[] vars) {
        try {
            // 更新距离
            for (int i = 0; i < vars.length; i++) {
                if (i % 2 == 0) {
                    stockyardList.get(i / 2).x = vars[i];
                } else {
                    stockyardList.get(i / 2).y = vars[i];
                }
            }
            // 计算距离矩阵
            double[][] distanceMatrix = new double[stockyardList.size()][constructionSiteList.size()];
            for (int i = 0; i < distanceMatrix.length; i++) {
                Stockyard stockyard = stockyardList.get(i);
                for (int j = 0; j < distanceMatrix[i].length; j++) {
                    ConstructionSite constructionSite = constructionSiteList.get(j);
                    distanceMatrix[i][j] = calcDistance(stockyard.x, stockyard.y, constructionSite.x, constructionSite.y);
                }
            }
            // 开始建模
//            cplex = new IloCplex();
            // 声明变量
            IloNumVar[][] x = new IloNumVar[stockyardList.size()][constructionSiteList.size()];
            for (int i = 0; i < x.length; i++) {
                for (int j = 0; j < x[i].length; j++) {
                    x[i][j] = cplex.numVar(0, Math.min(stockyardList.get(i).c, constructionSiteList.get(j).d));
                }
            }
            // 构造约束1：必须满足每个工地的需求
            for (int j = 0; j < constructionSiteList.size(); j++) {
                IloLinearNumExpr expr = cplex.linearNumExpr();
                for (int i = 0; i < x.length; i++) {
                    expr.addTerm(1, x[i][j]);
                }
                cplex.addEq(expr, constructionSiteList.get(j).d);
            }
            // 构造约束2：不能超过每个料场的储量
            for (int i = 0; i < stockyardList.size(); i++) {
                cplex.addLe(cplex.sum(x[i]), stockyardList.get(i).c);
            }
            // 声明目标函数
            IloLinearNumExpr target = cplex.linearNumExpr();
            for (int i = 0; i < x.length; i++) {
                for (int j = 0; j < x[i].length; j++) {
                    target.addTerm(distanceMatrix[i][j], x[i][j]);
                }
            }
            cplex.addMinimize(target);
            // 配置cplex
            cplex.setOut(null);
            cplex.setWarning(null);
            cplex.setParam(IloCplex.DoubleParam.EpOpt, EPS);
            // 开始求解
            double obj;
            if (cplex.solve()) {
                obj = cplex.getObjValue();
            } else {
                throw new RuntimeException("此题无解");
            }
            cplex.clearModel();
            return obj;
        } catch (Exception e) {
            throw new RuntimeException(e);
        }
    }

    // 复制蜘蛛猴
    SpiderMonkey copySpiderMonkey(SpiderMonkey old) {
        return new SpiderMonkey(old.curVars.clone(), old.curObjValue, old.fit);
    }

    public static void main(String[] args) throws IloException {
        // 料场列表
        List<Stockyard> stockyardList = new ArrayList<>();
        stockyardList.add(new Stockyard(0, 0, 20));
        stockyardList.add(new Stockyard(0, 0, 20));
        // 工地列表
        List<ConstructionSite> constructionSiteList = getInitConstructionSiteList();
        long s = System.currentTimeMillis();
        new SMO_Solver(stockyardList, constructionSiteList).solve();
        System.out.println("总用时: " + (System.currentTimeMillis() - s) / 1000d + " s");
    }

}

输出结果如下：

0 : 122.18394756340149
1 : 115.22137772493119
2 : 107.64015378143884
3 : 91.25562099863063
4 : 87.27034345689911
5 : 86.08239574989028
6 : 85.5990696901402
7 : 85.5675162336327
8 : 85.40967721594052
9 : 85.33079736839548
10 : 85.29533933219953
11 : 85.27770206804239
12 : 85.27194687104382
13 : 85.26919031512125
14 : 85.26862219676372
15 : 85.26750239234426
16 : 85.26663764799034
17 : 85.26636259651076
18 : 85.26628507350429
19 : 85.26618997863784
20 : 85.26605221867857
21 : 85.26605221867857
22 : 85.26604976682619
23 : 85.26604439739287
24 : 85.26604351700756
25 : 85.26604156579967
26 : 85.26604141563301
27 : 85.26604110466047
28 : 85.2660410153342
29 : 85.26604098426941
30 : 85.26604096946092
31 : 85.26604090826108
32 : 85.26604089012692
33 : 85.26604088253666
34 : 85.26604087662045
35 : 85.26604087515656
36 : 85.26604087344376
37 : 85.26604087292048
38 : 85.26604087263298
39 : 85.26604087239937
40 : 85.26604087223953
41 : 85.26604087214378
42 : 85.26604087212847
43 : 85.2660408721172
44 : 85.26604087211567
45 : 85.2660408721116
46 : 85.26604087211126
47 : 85.26604087211126
48 : 85.26604087211103
49 : 85.26604087211096
50 : 85.26604087211076
51 : 85.2660408721107
52 : 85.2660408721107
53 : 85.26604087211066
54 : 85.26604087211066
55 : 85.26604087211066
56 : 85.26604087211064
57 : 85.26604087211064
58 : 85.26604087211064
59 : 85.26604087211064
变量取值为：[7.25, 7.75, 3.2548825928531566, 5.652332445501583]
最优解为：85.26604087211064
总用时: 0.997 s