如何用Node.js的lp_solve求解带容量约束的工厂选址模型

问题背景 #

尝试复现带容量约束的工厂选址模型，Excel中目标函数通过SUMPRODUCT()实现，已用JavaScript编写sumProduct函数复现该逻辑，但使用lp_solve构建模型时，目标函数和容量约束出现错误，导致模型结果不符合预期。

复现的sumProduct函数 #

// Sample data arrays
/* const array1 = [1, 2, 3, 4, 5];
const array2 = [2, 3, 4, 5, 6]; */

// Function to calculate the sum product
export default function sumProduct(arr1, arr2) {
  if (arr1.length !== arr2.length) {
    throw new Error('Arrays must have the same length');
  }
  
  let result = 0;
  for (let i = 0; i < arr1.length; i++) {
    result += arr1[i] * arr2[i];
  }
  
  return result;
}

// Calculate the sum product
// const result = sumProduct(array1, array2); // OUTPUT: 70 (correct)

当前模型构建代码 #

import sumProduct  from "./sumProd.mjs";
import lp_solve from "lp_solve";

// Create a new LP model
const lp = new lp_solve.LinearProgram();

// Supply and demand locations
const locations = ["N_America", "S_America", "NCE", "MED", "APAC"];

// Transportation costs matrix
const transportationCosts = [
  [81, 92, 101, 130, 115],
  [117, 77, 108, 98, 100],
  [102, 105, 95, 119, 111],
  [115, 125, 90, 59, 74],
  [142, 100, 103, 105, 71],
];

// Fixed costs and capacities
const Fixed_Costs_L = [6000, 4500, 6500, 4100, 4000];
const Low_Capacity = [10, 10, 10, 10, 10];
const Fixed_Costs_H = [9000, 6750, 9750, 6150, 6000];
const High_Capacity = [20, 20, 20, 20, 20];

// Demand values for each location
const demand = [12, 8, 14, 16, 7];

// Create an object to store indices for demand locations
const demandIndices = {};
locations.forEach((location, index) => {
  demandIndices[location] = index;
});

// Add binary variables for opening/closing plants
const openPlants = {};
for (const location of locations) {
  openPlants[location] = lp.addColumn(`Open_${location}`, true, false);  // Indicate integer variable with 0/1 values
}


// Add columns (variables) for supply and demand allocations
const allocation = {};
for (const supplyLocation of locations) {
  allocation[supplyLocation] = {};

  for (const demandLocation of locations) {
    allocation[supplyLocation][demandLocation] = lp.addColumn(`${supplyLocation}_${demandLocation}`, true, true); // Indicate integer variable
  }

}

// Add constraints to ensure supply meets demand
for (const demandLocation of locations) {
  const demandConstraint = new lp_solve.Row();
  for (const supplyLocation of locations) {
    demandConstraint.Add(allocation[supplyLocation][demandLocation], 1);
  }
  lp.addConstraint(demandConstraint, "EQ", demand[demandIndices[demandLocation]], `Demand_${demandLocation}`);

}

// Add constraints for plant capacity based on the binary decision variable
for (const supplyLocation of locations) {
  const capacityConstraint = new lp_solve.Row();

  capacityConstraint.Add(openPlants[supplyLocation], High_Capacity[locations.indexOf(supplyLocation)]);
  capacityConstraint.Add(openPlants[supplyLocation], Low_Capacity[locations.indexOf(supplyLocation)]);
  
  for (const demandLocation of locations) {
    capacityConstraint.Subtract(allocation[supplyLocation][demandLocation], 1);
  }

  lp.addConstraint(capacityConstraint, "GE", 0, `Excess_Capacity_${supplyLocation}`);

}

// Add constraints to ensure supply and demand allocations are non-negative
for (const supplyLocation of locations) {
  for (const demandLocation of locations) {
    const nonNegativeConstraint = new lp_solve.Row();
    nonNegativeConstraint.Add(allocation[supplyLocation][demandLocation], 1);
    lp.addConstraint(nonNegativeConstraint, "GE", 0, `NonNegative_${supplyLocation}_${demandLocation}`);
  }
}

// Set the objective function to minimize total cost
const objective = new lp_solve.Row();
for (const supplyLocation of locations) {
  for (const demandLocation of locations) {
    objective.Add(allocation[supplyLocation][demandLocation], transportationCosts[locations.indexOf(supplyLocation)][locations.indexOf(demandLocation)]);
  }
  objective.Add(openPlants[supplyLocation], Fixed_Costs_L[locations.indexOf(supplyLocation)]);
  objective.Add(openPlants[supplyLocation], Fixed_Costs_H[locations.indexOf(supplyLocation)]);
}
lp.setObjective(objective, true); // true indicates minimizing

// Solve the LP optimization problem
const result = lp.solve();

console.log(lp.dumpProgram());
console.log(lp.solve());

// Display the results
if (result.description === "OPTIMAL") {
  console.log("Solution found:");
  for (const supplyLocation of locations) {
    for (const demandLocation of locations) {
      const value = lp.get(allocation[supplyLocation][demandLocation]);
      if (value > 0) {
        console.log(`${supplyLocation} -> ${demandLocation}: ${value}`);
      }
    }
    const openValue = lp.get(openPlants[supplyLocation]);
    console.log(`${supplyLocation} Open: ${openValue}`);
  }
  console.log(`Total Cost: $${lp.getObjectiveValue()}`);
} else {
  console.log("No feasible solution found.");
}

错误输出示例 #

目标函数（异常部分） #

minimize: +81 N_America_N_America +92 N_America_S_America +101 N_America_NCE +130 N_America_MED +115 N_America_APAC +15000 Open_N_America +117 S_America_N_America +77 S_America_S_America +108 S_America_NCE +98 S_America_MED +100 S_America_APAC +11250 Open_S_America +102 NCE_N_America +105 NCE_S_America +95 NCE_NCE +119 NCE_MED +111 NCE_APAC +16250 Open_NCE +115 MED_N_America +125 MED_S_America +90 MED_NCE +59 MED_MED +74 MED_APAC +10250 Open_MED +142 APAC_N_America +100 APAC_S_America +103 APAC_NCE +105 APAC_MED +71 APAC_APAC +10000 Open_APAC;

约束条件（容量约束异常部分） #

Excess_Capacity_N_America:  +30 Open_N_America -1 N_America_N_America -1 N_America_S_America -1 N_America_NCE -1 N_America_MED -1 N_America_APAC >= 0;
Excess_Capacity_S_America:  +30 Open_S_America -1 S_America_N_America -1 S_America_S_America -1 S_America_NCE -1 S_America_MED -1 S_America_APAC >= 0;
...

问题根源分析 #

1. 目标函数逻辑错误：将低容量和高容量的固定成本同时加到同一个Open_*变量上，导致固定成本系数为两者之和（比如Open_N_America的系数是6000+9000=15000），但实际每个工厂只能选择一种容量（低/高）或不开，不能同时承担两种固定成本。
2. 容量约束逻辑错误：将低容量和高容量的数值相加后乘以Open_*变量，得到30的系数，这不符合模型要求——工厂的容量应该是二选一，而非两者叠加。
3. 变量设计缺失：仅用一个二进制变量Open_*无法区分工厂是低容量运营还是高容量运营，需要拆分变量才能实现二选一的逻辑。

修正方案 #

1. 调整决策变量定义 #

为每个工厂创建两个二进制变量，分别对应低容量和高容量运营，并添加互斥约束：

// 替换原openPlants变量定义
const openPlantsL = {}; // 低容量运营
const openPlantsH = {}; // 高容量运营

for (const location of locations) {
  openPlantsL[location] = lp.addColumn(`Open_L_${location}`, true, false);
  openPlantsH[location] = lp.addColumn(`Open_H_${location}`, true, false);
  
  // 互斥约束：同一工厂不能同时开低容量和高容量
  const mutexConstraint = new lp_solve.Row();
  mutexConstraint.Add(openPlantsL[location], 1);
  mutexConstraint.Add(openPlantsH[location], 1);
  lp.addConstraint(mutexConstraint, "LE", 1, `Mutex_${location}`);
}

2. 修正容量约束 #

将容量约束改为低容量和高容量的线性组合，确保总分配量不超过所选容量：

// 替换原容量约束代码
for (const supplyLocation of locations) {
  const idx = locations.indexOf(supplyLocation);
  const capacityConstraint = new lp_solve.Row();
  
  // 低容量和高容量的贡献
  capacityConstraint.Add(openPlantsL[supplyLocation], Low_Capacity[idx]);
  capacityConstraint.Add(openPlantsH[supplyLocation], High_Capacity[idx]);
  
  // 减去总分配量
  for (const demandLocation of locations) {
    capacityConstraint.Subtract(allocation[supplyLocation][demandLocation], 1);
  }
  
  lp.addConstraint(capacityConstraint, "GE", 0, `Excess_Capacity_${supplyLocation}`);
}

3. 修正目标函数 #

固定成本部分改为低容量和高容量固定成本的线性组合：

// 替换原目标函数代码
const objective = new lp_solve.Row();
for (const supplyLocation of locations) {
  const idx = locations.indexOf(supplyLocation);
  
  // 运输成本部分不变
  for (const demandLocation of locations) {
    objective.Add(allocation[supplyLocation][demandLocation], transportationCosts[idx][locations.indexOf(demandLocation)]);
  }
  
  // 固定成本部分：低容量和高容量分别计算
  objective.Add(openPlantsL[supplyLocation], Fixed_Costs_L[idx]);
  objective.Add(openPlantsH[supplyLocation], Fixed_Costs_H[idx]);
}
lp.setObjective(objective, true);

4. 调整结果输出逻辑 #

输出时需要同时显示低容量和高容量的运营状态：

// 替换原结果输出代码
if (result.description === "OPTIMAL") {
  console.log("Solution found:");
  for (const supplyLocation of locations) {
    for (const demandLocation of locations) {
      const value = lp.get(allocation[supplyLocation][demandLocation]);
      if (value > 0) {
        console.log(`${supplyLocation} -> ${demandLocation}: ${value}`);
      }
    }
    const openLValue = lp.get(openPlantsL[supplyLocation]);
    const openHValue = lp.get(openPlantsH[supplyLocation]);
    console.log(`${supplyLocation} 低容量运营: ${openLValue}, 高容量运营: ${openHValue}`);
  }
  console.log(`Total Cost: $${lp.getObjectiveValue()}`);
} else {
  console.log("No feasible solution found.");
}

补充说明 #

关于sumProduct的疑问：lp_solve的目标函数和约束构建逻辑本质上就是线性组合，和SUMPRODUCT的计算逻辑一致，不需要额外调用自定义的sumProduct函数——通过Row对象的Add方法给变量赋予对应系数，最终模型会自动计算线性组合的结果，等价于Excel中SUMPRODUCT的作用。

内容的提问来源于stack exchange，提问作者Alex Vinulescu