问题
问题一求解
求解思路
该问题可以直接建立一个线性规划模型，然后使用cplex求解器来求解
模型
决策变量
代码实现
【存储结果的类】
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials; import lombok.Data; /** * @Author dam * @create 2023/8/15 20:32 */ @Data public class Solution { /** * 目标函数值 */ private double objectValue; /** * 料场i 向 料场j 运输的吨数 */ private double[][] xArr; }
【问题一模型求解主类】
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials; import ilog.concert.*; import ilog.cplex.IloCplex; /** * @Author dam * @create 2023/8/15 17:29 */ public class Question1 { public static void main(String[] args) throws IloException { 声明变量 // 各个工地的坐标位置 double[][] constructionSitePositionArr = { {1.25, 1.25}, {8.75, 0.75}, {0.5, 4.75}, {5.75, 5}, {3, 6.5}, {7.25, 7.75} }; // 各个工地的需求量 double[] demandArr = {3, 5, 4, 7, 6, 11}; // 各个料场的坐标位置 double[][] materialYardPositionArr = { {5, 1}, {2, 7} }; // 每个料场的水泥储量 double[] totalArr = {20, 20}; IloCplex cplex = new IloCplex(); Solution solution = cplexSolve(cplex, constructionSitePositionArr, demandArr, materialYardPositionArr, totalArr); // 问题求解完成，释放cplex资源 cplex.end(); System.out.println("目标函数：" + solution.getObjectValue()); for (int j = 0; j < solution.getXArr()[0].length; j++) { for (int i = 0; i < solution.getXArr().length; i++) { double value = solution.getXArr()[i][j]; if (value > 0) { System.out.println("料场" + (i + 1) + "分配给工地" + (j + 1) + "的吨数是：" + value); } } } } /** * cplex求解模型 * * @param constructionSitePositionArr * @param demandArr * @param materialYardPositionArr * @param totalArr * @return * @throws IloException */ public static Solution cplexSolve(IloCplex cplex, double[][] constructionSitePositionArr, double[] demandArr, double[][] materialYardPositionArr, double[] totalArr) throws IloException { // 存储求解结果 Solution solution = new Solution(); 建立模型 /// 决策变量 // 料场i 向 料场j 运输的吨数 IloNumVar[][] xArr = new IloNumVar[materialYardPositionArr.length][constructionSitePositionArr.length]; for (int i = 0; i < xArr.length; i++) { for (int j = 0; j < xArr[0].length; j++) { // 最大不超过需求量 xArr[i][j] = cplex.numVar(0, demandArr[j]); } } /// 添加目标函数：最小化运输费用，即最小化吨千米数 IloLinearNumExpr iloLinearNumExpr = cplex.linearNumExpr(); for (int i = 0; i < xArr.length; i++) { for (int j = 0; j < xArr[0].length; j++) { iloLinearNumExpr.addTerm(Question1.getDistance(materialYardPositionArr[i], constructionSitePositionArr[j]), xArr[i][j]); } } cplex.addMinimize(iloLinearNumExpr); /// 添加约束 // 约束1：每个工地收到的量等于自己的需求量 for (int j = 0; j < demandArr.length; j++) { IloLinearNumExpr iloNumExpr = cplex.linearNumExpr(); for (int i = 0; i < materialYardPositionArr.length; i++) { iloNumExpr.addTerm(xArr[i][j], 1); } cplex.addEq(iloNumExpr, demandArr[j]); } // 约束2：每个料场给工地所提供的水泥总量小于等于料场总量 for (int i = 0; i < totalArr.length; i++) { IloLinearNumExpr iloNumExpr = cplex.linearNumExpr(); for (int j = 0; j < constructionSitePositionArr.length; j++) { iloNumExpr.addTerm(xArr[i][j], 1); } cplex.addLe(iloNumExpr, totalArr[i]); } /// 模型求解 // 让cplex不要输出一些奇怪的东西，可以注释掉来看看效果 cplex.setOut(null); if (cplex.solve()) { solution.setObjectValue(cplex.getObjValue()); solution.setXArr(new double[xArr.length][xArr[0].length]); for (int j = 0; j < xArr[0].length; j++) { for (int i = 0; i < xArr.length; i++) { double value = cplex.getValue(xArr[i][j]); solution.getXArr()[i][j] = value; } } } // 清除模型，方便重复使用同一个cplex cplex.clearModel(); return solution; } /** * 获取两个位置之间的距离 * * @param position1 * @param position2 * @return */ private static double getDistance(double[] position1, double[] position2) { return Math.sqrt(Math.pow(position1[0] - position2[0], 2) + Math.pow(position1[1] - position2[1], 2)); } }
【运行结果】
目标函数：136.22751988318154 料场1分配给工地1的吨数是：3.0 料场1分配给工地2的吨数是：5.0 料场2分配给工地3的吨数是：4.0 料场1分配给工地4的吨数是：7.0 料场2分配给工地5的吨数是：6.0 料场1分配给工地6的吨数是：1.0 料场2分配给工地6的吨数是：10.0 Process finished with exit code 0
问题二求解
求解思路
如果直接使用建模的方式来求解，可以将两个料场的位置定为决策变量，即只需要修改上面的目标函数并增加决策变量，但是这样子建立的模型是非线性的，使用cplex无法求解，不过当然是可以使用其他工具来求解的，比如Matlab。
我这里是使用java语言来实现的，因此可以换一种方法来求解该问题，比如结合粒子群算法和上面的线性规划模型来求解。线性规划模型的作用是：给定料场的位置，可以求得料场给每个工地的供给量和总的目标函数值。而粒子群算法可以用来搜索料场的位置，将料场的位置定义为每个粒子的位置即可，在搜索的过程中，粒子每到一个新的位置，就可以调用模型来求解获得目标函数值，以此可以判断粒子的好坏。
代码实现
【克隆工具类】
因为在搜索的过程中，需要不断记录最优的料场位置，而本文的料场位置是使用二维数组来记录的（当然也可以使用一维数组，一维数组的代码实现更加方便，但是我个人觉得二维数组可以更好地理解），二维数组的clone方法不能直接使用，因此需要自己实现一个克隆的方法
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials.question2; /** * 克隆工具类 * @Author dam * @create 2023/8/16 10:10 */ public class CloneUtil { /** * 二维数组克隆 * @param oldArr * @return */ public static double[][] twoDArrClone(double[][] oldArr) { double[][] newArr = new double[oldArr.length][oldArr[0].length]; for (int i = 0; i < oldArr.length; i++) { newArr[i] = oldArr[i].clone(); } return newArr; } }
【粒子类】
每一个粒子需要保存一些个体信息，比如粒子当前所在位置、粒子速度、粒子曾经遍历到的历史最优解及其所在位置
初始化粒子的时候，需要设置粒子的位置，我这里直接让粒子随机在布局中的一个位置，这样可以让粒子群的粒子更加分散，找到质量较高的结果的概率更高
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials.question2; import lombok.Data; import java.util.Arrays; import java.util.Random; /** * 粒子类 * * @Author dam * @create 2023/8/15 20:40 */ @Data public class Particle { /** * 例子当前位置，即料场的位置 */ private double[][] materialYardPositionArr; /** * x轴方向上的速度 */ private double[] xVArr; /** * y轴方向上的速度 */ private double[] yVArr; /** * 该粒子找到的历史最优解 */ private double pBest; /** * 该粒子找到的历史最优解对应的位置 */ private double[][] bestMaterialYardPositionArr; public Particle(double xMin, double xMax, double yMin, double yMax, double xV, double yV, int materialYardNum, Random random) { this.xVArr = new double[materialYardNum]; Arrays.fill(this.xVArr, xV); this.yVArr = new double[materialYardNum]; Arrays.fill(this.yVArr, yV); // 初始化粒子信息，即初始化位置 this.materialYardPositionArr = new double[materialYardNum][2]; for (int i = 0; i < materialYardNum; i++) { this.materialYardPositionArr[i][0] = random.nextDouble() * (xMax - xMin) + xMin; this.materialYardPositionArr[i][1] = random.nextDouble() * (yMax - yMin) + yMin; } } /** * 设置最优解 */ public void setBest(double[][] materialYardPositionArr) { this.bestMaterialYardPositionArr = CloneUtil.twoDArrClone(materialYardPositionArr); } }
【粒子群算法类】
该类为粒子群算法类，主要用来初始化粒子群，不断迭代来控制粒子的速度变化及位置变化，在迭代过程中不断更新所找到的最优解
因为整个搜索过程需要调用多次cplex求解器来求解模型，我这里做了一点小优化，并不是每次求解模型都new一个IloCplex对象出来，而是只new一个，后面反复清空重新使用
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials.question2; import com.dam.intelligentWorkshopTrain.supply_of_raw_materials.Question1; import com.dam.intelligentWorkshopTrain.supply_of_raw_materials.Solution; import ilog.concert.IloException; import ilog.cplex.IloCplex; import java.util.Arrays; import java.util.Random; /** * 非线性模型，无法求解 * * @Author dam * @create 2023/8/15 19:59 */ public class PsoApi { /** * 粒子数量 */ private int particleNum; /** * 个体学习因子，设置得越大，粒子越容易根据自己的想法飞行，若设置过大，容易跳出局部最优，但收敛较慢 */ private double c1; /** * 社会学习因子，设置得越大，粒子越容易根据群体的想法飞行，若设置过大，容易陷入局部最优，收敛较快 */ private double c2; /** * 速度最大值 */ private double vMax; /** * 速度的惯性权重 */ private double w; /** * 迭代次数 */ private int genMax; public PsoApi(int particleNum, double c1, double c2, double vMax, double w, int genMax) { this.particleNum = particleNum; this.c1 = c1; this.c2 = c2; this.vMax = vMax; this.w = w; this.genMax = genMax; } /** * 求解 */ public void solve(double[][] constructionSitePositionArr, double[] demandArr, double[] totalArr) throws IloException { 变量声明 // 存储粒子 Particle[] particleArr; // 所有粒子找到的最优解（由于问题为最小化问题，设置初始最优值为较大的数） double globalBest = Double.MAX_VALUE; // 群体最优解对应的x和y double[][] bestMaterialYardPositionArr = null; // 随机数工具 Random random = new Random(); // cplex对象 IloCplex cplex = new IloCplex(); long start = System.currentTimeMillis(); int materialYardNum = totalArr.length; 初始化粒子 // 根据工地找出最大最小坐标 double xMin = Double.MAX_VALUE; double xMax = -Double.MAX_VALUE; double yMin = Double.MAX_VALUE; double yMax = -Double.MAX_VALUE; for (int i = 0; i < constructionSitePositionArr.length; i++) { double curX = constructionSitePositionArr[i][0]; double curY = constructionSitePositionArr[i][1]; xMin = Math.min(curX, xMin); xMax = Math.max(curX, xMax); yMin = Math.min(curY, yMin); yMax = Math.max(curY, yMax); } particleArr = new Particle[this.particleNum]; for (int i = 0; i < particleArr.length; i++) { // 初始化粒子群，注意：这里设置每个粒子的速度一样，读者可以根据自己的喜爱进行设置 // 初始化粒子的坐标和速度 particleArr[i] = new Particle(xMin, xMax, yMin, yMax, 0.01, 0.01, materialYardNum, random); // 初始化粒子的函数值(粒子还没有开始飞，当前位置肯定是找到过的最优位置啦) particleArr[i].setBest(particleArr[i].getMaterialYardPositionArr()); Solution solution = Question1.cplexSolve(cplex, constructionSitePositionArr, demandArr, particleArr[i].getMaterialYardPositionArr(), totalArr); particleArr[i].setPBest(solution.getObjectValue()); // 由于问题为最小化问题，目标函数越小越好 if (solution.getObjectValue() < globalBest) { bestMaterialYardPositionArr = CloneUtil.twoDArrClone(particleArr[i].getMaterialYardPositionArr()); globalBest = solution.getObjectValue(); } } 开始求解 for (int i = 0; i < this.genMax; i++) { // 对每个粒子进行操作 for (int j = 0; j < this.particleNum; j++) { // 更新速度 updateSpeed(particleArr, bestMaterialYardPositionArr, random, materialYardNum, j); // 更新位置 updatePosition(materialYardNum, xMin, xMax, yMin, yMax, particleArr[j]); /// 更新粒子历史最优解和粒子全体最优解 Solution solution = Question1.cplexSolve(cplex, constructionSitePositionArr, demandArr, particleArr[j].getMaterialYardPositionArr(), totalArr); if (solution.getObjectValue() < particleArr[j].getPBest()) { particleArr[j].setBest(particleArr[j].getMaterialYardPositionArr()); particleArr[j].setPBest(solution.getObjectValue()); } // 由于问题为最小化问题，目标函数越小越好 if (solution.getObjectValue() < globalBest) { bestMaterialYardPositionArr = CloneUtil.twoDArrClone(particleArr[j].getMaterialYardPositionArr()); globalBest = solution.getObjectValue(); System.out.println("当前最优解：" + globalBest); System.out.println("料场最优坐标"); for (int w = 0; w < bestMaterialYardPositionArr.length; w++) { System.out.println(Arrays.toString(bestMaterialYardPositionArr[w])); } System.out.println("当前是第" + (i + 1) + "代"); System.out.println(); } } } // 问题求解完成，释放cplex资源 cplex.end(); // 输出保留6位小数 System.out.println("---------------------------结果输出---------------------------"); System.out.println("最优目标函数值：" + String.format("%.6f", globalBest)); System.out.println("料场最优坐标"); for (int i = 0; i < bestMaterialYardPositionArr.length; i++) { System.out.println(Arrays.toString(bestMaterialYardPositionArr[i])); } System.out.println("求解时间：" + (System.currentTimeMillis() - start) + "ms"); } /** * 更新例子的位置 * @param materialYardNum * @param xMin * @param xMax * @param yMin * @param yMax * @param particleArr */ private void updatePosition(int materialYardNum, double xMin, double xMax, double yMin, double yMax, Particle particleArr) { for (int k = 0; k < materialYardNum; k++) { //todo 这部分的冗余代码比较多，读者可以想办法优化，比如将xVArr和yVArr合并为一个二维的数组，这样一个遍历就可以解决了 double nextX = particleArr.getMaterialYardPositionArr()[k][0] + particleArr.getXVArr()[k]; // 处理越界 if (nextX > xMax) { nextX = xMax; } else if (nextX < xMin) { nextX = xMin; } particleArr.getMaterialYardPositionArr()[k][0] = nextX; double nextY = particleArr.getMaterialYardPositionArr()[k][1] + particleArr.getYVArr()[k]; // 处理越界 if (nextY > yMax) { nextY = yMax; } else if (nextY < yMin) { nextY = yMin; } particleArr.getMaterialYardPositionArr()[k][1] = nextY; } } /** * 更新粒子的速度 * @param particleArr * @param bestMaterialYardPositionArr * @param random * @param materialYardNum * @param j */ private void updateSpeed(Particle[] particleArr, double[][] bestMaterialYardPositionArr, Random random, int materialYardNum, int j) { for (int k = 0; k < materialYardNum; k++) { //todo 这部分的冗余代码比较多，读者可以想办法优化，比如将xVArr和yVArr合并为一个二维的数组，这样一个遍历就可以解决了 // 更新x轴方向上的速度 particleArr[j].getXVArr()[k] = this.w * particleArr[j].getXVArr()[k] + this.c1 * random.nextDouble() * (particleArr[j].getBestMaterialYardPositionArr()[k][0] - particleArr[j].getMaterialYardPositionArr()[k][0]) + this.c2 * random.nextDouble() * (bestMaterialYardPositionArr[k][0] - particleArr[j].getMaterialYardPositionArr()[k][0]); // 处理x轴方向速度越界 if (particleArr[j].getXVArr()[k] > this.vMax) { particleArr[j].getXVArr()[k] = this.vMax; } else if (particleArr[j].getXVArr()[k] < -this.vMax) { particleArr[j].getXVArr()[k] = -this.vMax; } // 更新y轴方向上的速度 particleArr[j].getYVArr()[k] = this.w * particleArr[j].getYVArr()[k] + this.c1 * random.nextDouble() * (particleArr[j].getBestMaterialYardPositionArr()[k][1] - particleArr[j].getMaterialYardPositionArr()[k][1]) + this.c2 * random.nextDouble() * (bestMaterialYardPositionArr[k][1] - particleArr[j].getMaterialYardPositionArr()[k][1]); // 处理y轴方向速度越界 if (particleArr[j].getYVArr()[k] > this.vMax) { particleArr[j].getYVArr()[k] = this.vMax; } else if (particleArr[j].getYVArr()[k] < -this.vMax) { particleArr[j].getYVArr()[k] = -this.vMax; } } } }
【主类】
该类主要用来声明数据，给粒子群算法设置参数，最后调用粒子群算法来对问题进行求解
package com.dam.intelligentWorkshopTrain.supply_of_raw_materials.question2; import ilog.concert.IloException; /** * @Author dam * @create 2023/8/15 21:24 */ public class Question2_PSO { public static void main(String[] args) { 声明变量 // 各个工地的坐标位置 double[][] constructionSitePositionArr = { {1.25, 1.25}, {8.75, 0.75}, {0.5, 4.75}, {5.75, 5}, {3, 6.5}, {7.25, 7.75} }; // 各个工地的需求量 double[] demandArr = {3, 5, 4, 7, 6, 11}; // 每个料场的水泥储量 double[] totalArr = {20, 20}; // 设置粒子群算法的参数 PsoApi psoApi = new PsoApi(100, 0.05, 0.05, 5, 0.9, 100); try { // 调用粒子群算法 psoApi.solve(constructionSitePositionArr, demandArr, totalArr); } catch (IloException e) { throw new RuntimeException(e); } } }
【运行】
当前最优解：88.50017762115978 料场最优坐标 [6.995460835765463, 7.431965294573079] [2.4232613334428765, 5.420225506839392] 当前是第1代 当前最优解：88.43200001652829 料场最优坐标 [7.003560835765462, 7.4400652945730785] [2.4313613334428767, 5.428325506839392] 当前是第2代 当前最优解：88.37092613583344 料场最优坐标 [7.010850835765463, 7.447355294573079] [2.4386513334428765, 5.435615506839392] 当前是第3代 当前最优解：88.31619438714196 料场最优坐标 [7.017411835765462, 7.453916294573078] [2.4452123334428766, 5.442176506839392] 当前是第4代 当前最优解：88.26712846008976 料场最优坐标 [7.023316735765462, 7.459821194573078] [2.4511172334428766, 5.448081406839392] 当前是第5代 当前最优解：88.22312716663495 料场最优坐标 [7.028631145765463, 7.465135604573079] [2.4564316434428766, 5.453395816839392] 当前是第6代 当前最优解：88.1836555959095 料场最优坐标 [7.0334141147654625, 7.469918573573079] [2.4612146124428764, 5.458178785839392] 当前是第7代 当前最优解：88.14823739455572 料场最优坐标 [7.0377187868654625, 7.474223245673079] [2.4655192845428764, 5.462483457939392] 当前是第8代 当前最优解：87.87110543103842 料场最优坐标 [7.37221941415926, 7.472158486930282] [3.744400991411965, 5.193926703167144] 当前是第9代 当前最优解：86.94994195734931 料场最优坐标 [7.1143623064881805, 7.627892453016968] [3.4019200051871716, 5.430391345444132] 当前是第10代 当前最优解：86.23975720549879 料场最优坐标 [7.324184226794703, 7.663656711228883] [3.012995154097673, 5.642541692214598] 当前是第11代 当前最优解：86.23745998092613 料场最优坐标 [7.236137326506472, 7.75] [2.804696509875521, 5.792557866752227] 当前是第13代 当前最优解：85.79420546873902 料场最优坐标 [7.277291908434675, 7.75] [3.042386487841036, 5.490155588705559] 当前是第13代 当前最优解：85.7401223802724 料场最优坐标 [7.252795182230648, 7.75] [2.9139246554783913, 5.622177963512921] 当前是第15代 当前最优解：85.71122420252955 料场最优坐标 [7.2226019632563165, 7.75] [3.2656022819683046, 5.871079972622107] 当前是第19代 当前最优解：85.49458063165437 料场最优坐标 [7.245769701653424, 7.75] [3.2946251539561096, 5.904325759516212] 当前是第20代 当前最优解：85.4509241798177 料场最优坐标 [7.264671022225204, 7.75] [3.3014973294756103, 5.749172090266563] 当前是第25代 当前最优解：85.383103191951 料场最优坐标 [7.258010737754079, 7.75] [3.340243912257952, 5.695326263970876] 当前是第25代 当前最优解：85.37611196230958 料场最优坐标 [7.260859212787864, 7.75] [3.2305121221006603, 5.668446593189184] 当前是第28代 当前最优解：85.37028591165853 料场最优坐标 [7.2536923453472015, 7.75] [3.376598914930889, 5.6973535637469315] 当前是第29代 当前最优解：85.32396377315396 料场最优坐标 [7.252168191046689, 7.75] [3.1676286576076826, 5.760573217144864] 当前是第29代 当前最优解：85.30131868341802 料场最优坐标 [7.249047625903025, 7.75] [3.20787665106796, 5.760005370791448] 当前是第31代 当前最优解：85.27628292994825 料场最优坐标 [7.250787629818613, 7.75] [3.263808468189787, 5.679476503001795] 当前是第41代 当前最优解：85.27025155292363 料场最优坐标 [7.249905732152565, 7.75] [3.268868570723304, 5.679396795227463] 当前是第45代 当前最优解：85.26812668241098 料场最优坐标 [7.250191427294184, 7.75] [3.2498860087141708, 5.661247085205052] 当前是第57代 当前最优解：85.26657211227986 料场最优坐标 [7.249971762881258, 7.75] [3.24862404799676, 5.660537597317892] 当前是第71代 当前最优解：85.26646536552198 料场最优坐标 [7.250041161365652, 7.75] [3.252702374447127, 5.653289733955173] 当前是第90代 当前最优解：85.2662612100008 料场最优坐标 [7.250020608545754, 7.75] [3.2528207689360498, 5.654146585813984] 当前是第91代 当前最优解：85.26608067383026 料场最优坐标 [7.250002111007846, 7.75] [3.2529273239760803, 5.654917752486913] 当前是第92代 ---------------------------结果输出--------------------------- 最优目标函数值：85.266081 料场最优坐标 [7.250002111007846, 7.75] [3.2529273239760803, 5.654917752486913] 求解时间：956ms Process finished with exit code 0
