Pooling
Problem
In the Pooling Problem, flows of raw materials with specific attribute concentrations are blended to create final products. The attribute concentrations of the final products must be within tolerance intervals. Indeed, these attributes are interesting primary constituents of the materials. The flows can either be sent directly from the source nodes to the target nodes or blended at intermediate pools and then sent to the target nodes. In addition, the capacities of the flow graph and of the different nodes must be respected, and the raw materials supply cannot be exceeded. The objective is to maximize the operation’s profit. The gains come from selling the final products, while the costs correspond to the raw materials to buy.
In other words, the Pooling Problem combines aspects of the Minimal Cost Flow Problem and the Blending Problem.
Principles learned
- Use the JSON module to read the input files
- Add float decision variables to model the flows
- Use non-linear operators
Data
The instances come from this GitHub repository and follow the JSON format:
- “components”: array of objects characterizing the raw materials with:
- “name”: name of the source node
- “lower”: minimal outflow (unused as equal to 0)
- “upper”: maximal outflow (supply)
- “price”: price for one unit of the raw material
- “quality”: object containing the concentration of each attribute in the raw material
- “products”: array of objects characterizing the final products with:
- “name”: name of the target node
- “lower”: minimal inflow (demand)
- “upper”: maximal inflow (capacity)
- “price”: sale price for one unit of the final product
- “quality_lower”: object containing the minimum concentration of each attribute in the final product
- “quality_upper”: object containing the maximum concentration of each attribute in the final product
- “pool_size”: object containing the capacities of the pools
- “component_to_product_bound”: array of objects characterizing the edges between the source nodes and the target nodes with:
- “component”: name of the source node c
- “product”: name of the target node p
- “bound”: maximum flow on edge (c, p)
- “cost”: cost of one unit of flow on edge (c, p)
- “component_to_pool_fraction”: array of objects characterizing the edges between the source node and the pool nodes with:
- “component”: name of the source node c
- “pool”: name of the pool node o
- “fraction”: maximum proportion of inflow at pool p coming from component c
- “cost”: cost of one unit of flow in edge (c, o)
- “pool_to_product_bound”: array of objects characterizing the edges between the pool nodes and the target nodes with:
- “pool”: name of the pool node o
- “product”: name of the target node p
- “bound”: maximum flow on edge (o, p)
- “cost”: cost of one unit of flow in edge (o, p)
The price of the raw materials and of the final products is either conveyed in the “price” or “cost” feature on the edges. The “cost” representation is used for the “randstd” instances, while the “price” representation is used for the others.
We use JSON libraries for each API model to read the instance data and write the output solution. For example, we use “Newtonsoft.Json” for C#, “gson-2.8.8” for Java, “json” for Python, and “nlohmann/json.hpp” for C++. For Hexaly Modeler, we use the JSON module.
Program
The Hexaly model for the Pooling Problem corresponds to the Q formulation of the problem, which uses both flow proportions and flow quantities. We start by declaring three arrays of float decision variables, representing:
- The flow from each source node c to each target node p;
- The flow from each pool o to each target node p;
- The proportion of the inflow at each pool o coming from each source node c.
By multiplying the two latter quantities, we compute the flow coming from each source node c to each target node p through each pool o.
We can then start writing the constraints. Firstly, we ensure that the sum of the proportions of inflow at each pool from the source nodes is equal to 1. Then, the total inflow at each target node must satisfy the demand while respecting the product capacity. Similarly, the total outflow at each source node cannot exceed the supply of the raw materials. Finally, the concentration of each attribute for each final product must be within a certain tolerance interval. To enforce this constraint, we compute the attribute quantities coming directly from the source nodes and from the pools, as well as the total inflow at the target node.
Lastly, we compute the value of the objective function. Total profit is equal to the total revenue acquired by selling the products minus the cost of the initial components minus the flow costs.
- Execution
-
hexaly pooling.hxm inFileName=instances/haverly1.json [hxTimeLimit=] [solFileName=]
use json;
/* Read instance data */
function input() {
local usage = "Usage: hexaly pooling.hxm "
+ "inFileName=inputFile [hxTimeLimit=timeLimit]";
if (inFileName == nil) throw usage;
problem = json.parse(inFileName);
nbComponents = problem.components.count();
nbProducts = problem.products.count();
nbPools = problem.pool_size.count();
nbAttributes = problem.components[0].quality.count();
// Components
componentQualities[c in 0...nbComponents] = problem.components[c].quality.values();
for [c in 0...nbComponents]
componentIdx[problem.components[c].name] = c;
// Final products (blendings)
for [p in 0...nbProducts]
productIdx[problem.products[p].name] = p;
for [p in 0...nbProducts] {
if (problem.products[p]["quality_lower"] == nil)
minTolerance[p][k in 0...nbAttributes] = 0;
else
minTolerance[p] = problem.products[p].quality_lower.values();
}
maxTolerance[p in 0...nbProducts] = problem.products[p].quality_upper.values();
// Intermediate pools
poolCapacities = problem.pool_size.values();
poolNames = problem.pool_size.keys();
for [o in 0...nbPools]
poolIdx[poolNames[o]] = o;
// Flow graph
// Edges from the components to the products
upperBoundComponentToProduct[c in 0...nbComponents][p in 0...nbProducts] = 0;
costComponentToProduct[c in 0...nbComponents][p in 0...nbProducts] = 0;
// Edges from the components to the pools
upperBoundFractionComponentToPool[c in 0...nbComponents][o in 0...nbPools] = 0;
costComponentToPool[c in 0...nbComponents][o in 0...nbPools] = 0;
// Edges from the pools to the products
upperBoundPoolToProduct[o in 0...nbPools][p in 0...nbProducts] = 0;
costPoolToProduct[o in 0...nbPools][p in 0...nbProducts] = 0;
// Bound and cost on the edges
for [edge in problem.component_to_product_bound] {
upperBoundComponentToProduct[componentIdx[edge.component]]
[productIdx[edge.product]] = edge.bound;
if (count(edge) > 3)
costComponentToProduct[componentIdx[edge.component]]
[productIdx[edge.product]] = edge.cost;
}
for [edge in problem.component_to_pool_fraction] {
upperBoundFractionComponentToPool[componentIdx[edge.component]]
[poolIdx[edge.pool]] = edge.fraction;
if (count(edge) > 3)
costComponentToPool[componentIdx[edge.component]]
[poolIdx[edge.pool]] = edge.cost;
}
for [edge in problem.pool_to_product_bound] {
upperBoundPoolToProduct[poolIdx[edge.pool]]
[productIdx[edge.product]] = edge.bound;
if (count(edge) > 3)
costPoolToProduct[poolIdx[edge.pool]]
[productIdx[edge.product]] = edge.cost;
}
}
/* Declare the optimization model */
function model() {
// Decision variables
// Flow from the components to the products
flowComponentToProduct[c in 0...nbComponents][p in 0...nbProducts] <-
float(0, upperBoundComponentToProduct[c][p]);
// Fraction of the total flow in pool o coming from the component c
fractionComponentToPool[c in 0...nbComponents][o in 0...nbPools] <-
float(0, upperBoundFractionComponentToPool[c][o]);
// Flow from the pools to the products
flowPoolToProduct[o in 0...nbPools][p in 0...nbProducts] <-
float(0, upperBoundPoolToProduct[o][p]);
// Flow from the components to the products and passing by the pools
flowComponentToProductByPool[c in 0...nbComponents][o in 0...nbPools][p in 0...nbProducts] <-
fractionComponentToPool[c][o] * flowPoolToProduct[o][p];
// Constraints
// Proportion
for [o in 0...nbPools] {
proportion <- sum[c in 0...nbComponents](fractionComponentToPool[c][o]);
constraint proportion == 1;
}
// Component supply
for [c in 0...nbComponents] {
flowToProducts <- sum[p in 0...nbProducts](flowComponentToProduct[c][p]);
flowToPools <- sum[o in 0...nbPools][p in 0...nbProducts]
(flowComponentToProductByPool[c][o][p]);
totalOutFlow <- flowToPools + flowToProducts;
constraint totalOutFlow <= problem.components[c].upper;
}
// Pool capacity (bounds on edges)
for [c in 0...nbComponents] {
for [o in 0...nbPools] {
flowComponentToPool <- sum[p in 0...nbProducts]
(flowComponentToProductByPool[c][o][p]);
edgeCapacity <- poolCapacities[o] * fractionComponentToPool[c][o];
constraint flowComponentToPool <= edgeCapacity;
}
}
// Product capacity
for [p in 0...nbProducts] {
flowFromPools <- sum[o in 0...nbPools](flowPoolToProduct[o][p]);
flowFromComponents <- sum[c in 0...nbComponents](flowComponentToProduct[c][p]);
totalInFlow <- flowFromComponents + flowFromPools;
constraint totalInFlow <= problem.products[p].upper;
constraint totalInFlow >= problem.products[p].lower;
}
// Product tolerance
for [p in 0...nbProducts] {
for [k in 0...nbAttributes] {
// Attribute from the components
attributeFromComponents <- sum[c in 0...nbComponents]
(componentQualities[c][k] * flowComponentToProduct[c][p]);
// Attribute from the pools
attributeFromPools <- sum[c in 0...nbComponents](componentQualities[c][k] *
sum[o in 0...nbPools](flowComponentToProductByPool[c][o][p]));
// Total flow in the blending
totalFlowIn <- sum[o in 0...nbPools](flowPoolToProduct[o][p]) +
sum[c in 0...nbComponents](flowComponentToProduct[c][p]);
totalAttributeIn <- attributeFromComponents + attributeFromPools;
constraint totalAttributeIn >= minTolerance[p][k] * totalFlowIn;
constraint totalAttributeIn <= maxTolerance[p][k] * totalFlowIn;
}
}
// Objective function
// Cost of the flows from the components directly to the products
directFlowCost <- sum[c in 0...nbComponents][p in 0...nbProducts]
(flowComponentToProduct[c][p] * costComponentToProduct[c][p]);
// Cost of the flows from the components to the products passing by the pools
undirectFlowCost <- sum[c in 0...nbComponents][o in 0...nbPools][p in 0...nbProducts]
((costComponentToPool[c][o] + costPoolToProduct[o][p]) *
flowComponentToProductByPool[c][o][p]);
// Gain of selling the final products
productsGain <- sum[p in 0...nbProducts](problem.products[p].price *
(sum[o in 0...nbPools](flowPoolToProduct[o][p]) +
sum[c in 0...nbComponents](flowComponentToProduct[c][p])));
// Cost of buying the components
componentsCost <- sum[c in 0...nbComponents](problem.components[c].price *
(sum[o in 0...nbPools][p in 0...nbProducts]
(fractionComponentToPool[c][o] * flowPoolToProduct[o][p]) +
sum[p in 0...nbProducts](flowComponentToProduct[c][p])));
profit <- productsGain - componentsCost - (directFlowCost + undirectFlowCost);
// Maximize the total profit
maximize profit;
}
/* Parameterize the solver */
function param() {
if (hxTimeLimit == nil) hxTimeLimit = 20;
}
/* Write the solution in a file */
function output() {
if (solFileName == nil) return;
// Solution flow from the components to the products
for [c in 0...nbComponents] {
for [p in 0...nbProducts] {
component_to_poduct[c * nbProducts + p] = {
"component" = problem.components[c].name,
"product" = problem.products[p].name,
"flow" = flowComponentToProduct[c][p].value};
}
}
// Solution fraction of the inflow at pool o coming from the component c
for [c in 0...nbComponents] {
for [o in 0...nbPools] {
component_to_pool_fraction[c * nbPools + o] = {
"component" = problem.components[c].name,
"pool" = poolNames[o],
"flow" = fractionComponentToPool[c][o].value};
}
}
// Solution flows from the pools to the products
for [o in 0...nbPools] {
for [p in 0...nbProducts] {
pool_to_product[o * nbProducts + p] = {
"pool" = poolNames[o],
"product" = problem.products[p].name,
"flow" = flowPoolToProduct[o][p].value};
}
}
json.dump({"objective" = profit.value, "solution" = {
"component_to_pool_fraction" = component_to_pool_fraction,
"component_to_product" = component_to_poduct ,
"pool_to_product" = pool_to_product}}, solFileName);
}
- Execution (Windows)
-
set PYTHONPATH=%HX_HOME%\bin\pythonpython pooling.py instances\haverly1.json
- Execution (Linux)
-
export PYTHONPATH=/opt/hexaly_13_0/bin/pythonpython pooling.py instances\haverly1.json
import hexaly.optimizer
import json
import sys
class PoolingInstance:
#
# Read instance data
#
def __init__(self, instance_file):
with open(instance_file) as problem:
problem = json.load(problem)
self.nbComponents = len(problem["components"])
self.nbProducts = len(problem["products"])
self.nbAttributes = len(problem["components"][0]["quality"])
self.nbPools = len(problem["pool_size"])
# Components
self.componentPrices = [problem["components"][c]["price"]
for c in range(self.nbComponents)]
self.componentSupplies = [problem["components"][c]["upper"]
for c in range(self.nbComponents)]
self.componentQuality = [list(problem["components"][c]["quality"].values())
for c in range(self.nbComponents)]
self.componentNames = [problem["components"][c]["name"]
for c in range(self.nbComponents)]
componentsIdx = {}
for c in range(self.nbComponents):
componentsIdx[problem["components"][c]["name"]] = c
# Final products (blendings)
self.productPrices = [problem["products"][p]["price"]
for p in range(self.nbProducts)]
self.productCapacities = [problem["products"][p]["upper"]
for p in range(self.nbProducts)]
self.demand = [problem["products"][p]["lower"]
for p in range(self.nbProducts)]
self.productNames = [problem["products"][p]["name"]
for p in range(self.nbProducts)]
productIdx = {}
for p in range(self.nbProducts):
productIdx[problem["products"][p]["name"]] = p
self.minTolerance = [[0 for _ in range(self.nbAttributes)]
if (problem["products"][p]["quality_lower"] == None)
else list(problem["products"][p]["quality_lower"].values())
for p in range(self.nbProducts)]
self.maxTolerance = [list(problem["products"][p]["quality_upper"].values())
for p in range(self.nbProducts)]
# Intermediate pools
self.poolNames = list(problem["pool_size"].keys())
self.poolCapacities = [problem["pool_size"][o] for o in self.poolNames]
poolIdx = {}
for o in range(self.nbPools):
poolIdx[self.poolNames[o]] = o
# Flow graph
# Edges from the components to the products
self.upperBoundComponentToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbComponents)]
self.costComponentToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbComponents)]
# Edges from the components to the pools
self.upperBoundFractionComponentToPool = [[0 for _ in range(self.nbPools)]
for _ in range(self.nbComponents)]
self.costComponentToPool = [[0 for _ in range(self.nbPools)]
for _ in range(self.nbComponents)]
# Edges from the pools to the products
self.upperBoundPoolToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbPools)]
self.costPoolToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbPools)]
# Bound and cost on the edges
for edge in problem["component_to_product_bound"]:
self.upperBoundComponentToProduct[componentsIdx[edge["component"]]] \
[productIdx[edge["product"]]] = edge["bound"]
if len(edge) > 3:
self.costComponentToProduct[componentsIdx[edge["component"]]] \
[productIdx[edge["product"]]] = edge["cost"]
for edge in problem["component_to_pool_fraction"]:
self.upperBoundFractionComponentToPool[componentsIdx[edge["component"]]] \
[poolIdx[edge["pool"]]] = edge["fraction"]
if len(edge) > 3:
self.costComponentToPool[componentsIdx[edge["component"]]] \
[poolIdx[edge["pool"]]] = edge["cost"]
for edge in problem["pool_to_product_bound"]:
self.upperBoundPoolToProduct[poolIdx[edge["pool"]]] \
[productIdx[edge["product"]]] = edge["bound"]
if len(edge) > 3:
self.costPoolToProduct[poolIdx[edge["pool"]]] \
[productIdx[edge["product"]]] = edge["cost"]
def main(instance_file, output_file, time_limit):
data = PoolingInstance(instance_file)
with hexaly.optimizer.HexalyOptimizer() as optimizer:
#
# Declare the optimization model
#
model = optimizer.model
# Decision variables
# Flow from the components to the products
flowComponentToProduct = [[model.float(0, data.upperBoundComponentToProduct[c][p])
for p in range(data.nbProducts)] for c in range(data.nbComponents)]
# Fraction of the total flow in pool o coming from the component c
fractionComponentToPool = [[model.float(0, data.upperBoundFractionComponentToPool[c][o])
for o in range(data.nbPools)]
for c in range(data.nbComponents)]
# Flow from the pools to the products
flowPoolToProduct = [
[model.float(0, data.upperBoundPoolToProduct[o][p])
for p in range(data.nbProducts)] for o in range(data.nbPools)]
# Flow from the components to the products and passing by the pools
flowComponentToProductByPool = [
[[fractionComponentToPool[c][o] * flowPoolToProduct[o][p]
for p in range(data.nbProducts)] for o in range(data.nbPools)]
for c in range(data.nbComponents)]
# Constraints
# Proportion
for o in range(data.nbPools):
proportion = model.sum(fractionComponentToPool[c][o]
for c in range(data.nbComponents))
model.constraint(proportion == 1)
# Component supply
for c in range(data.nbComponents):
flowToProducts = model.sum(flowComponentToProduct[c][p]
for p in range(data.nbProducts))
flowToPools = model.sum(flowComponentToProductByPool[c][o][p]
for p in range(data.nbProducts) for o in range(data.nbPools))
totalOutFlow = model.sum(flowToPools, flowToProducts)
model.constraint(totalOutFlow <= data.componentSupplies[c])
# Pool capacity (bounds on edges)
for c in range(data.nbComponents):
for o in range(data.nbPools):
flowComponentToPool = model.sum(flowComponentToProductByPool[c][o][p]
for p in range(data.nbProducts))
edgeCapacity = model.prod(data.poolCapacities[o],
fractionComponentToPool[c][o])
model.constraint(flowComponentToPool <= edgeCapacity)
# Product capacity
for p in range(data.nbProducts):
flowFromPools = model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))
flowFromComponents = model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents))
totalInFlow = model.sum(flowFromComponents, flowFromPools)
model.constraint(totalInFlow <= data.productCapacities[p])
model.constraint(totalInFlow >= data.demand[p])
# Product tolerance
for p in range(data.nbProducts):
for k in range(data.nbAttributes):
# Attribute from the components
attributeFromComponents = model.sum(
data.componentQuality[c][k] * flowComponentToProduct[c][p]
for c in range(data.nbComponents))
# Attribute from the pools
attributeFromPools = model.sum(
data.componentQuality[c][k] * flowComponentToProductByPool[c][o][p]
for o in range(data.nbPools) for c in range(data.nbComponents))
# Total flow in the blending
totalFlowIn = model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents)) \
+ model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))
totalAttributeIn = model.sum(attributeFromComponents, attributeFromPools)
model.constraint(totalAttributeIn >= data.minTolerance[p][k] * totalFlowIn)
model.constraint(totalAttributeIn <= data.maxTolerance[p][k] * totalFlowIn)
# Objective function
# Cost of the flows from the components directly to the products
directFlowCost = model.sum(
data.costComponentToProduct[c][p] * flowComponentToProduct[c][p]
for c in range(data.nbComponents) for p in range(data.nbProducts))
# Cost of the flows from the components to the products passing by the pools
undirectFlowCost = model.sum(
(data.costComponentToPool[c][o] + data.costPoolToProduct[o][p]) *
flowComponentToProductByPool[c][o][p]
for c in range(data.nbComponents)
for o in range(data.nbPools) for p in range(data.nbProducts))
# Gain of selling the final products
productsGain = model.sum((model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents))
+ model.sum(flowPoolToProduct[o][p]
for o in range(data.nbPools)))
* data.productPrices[p] for p in range(data.nbProducts))
# Cost of buying the components
componentsCost = model.sum(
(model.sum(flowComponentToProduct[c][p]
for p in range(data.nbProducts))
+ model.sum(fractionComponentToPool[c][o] * flowPoolToProduct[o][p]
for p in range(data.nbProducts)
for o in range(data.nbPools)))
* data.componentPrices[c] for c in range(data.nbComponents))
profit = productsGain - componentsCost - (directFlowCost + undirectFlowCost)
# Maximize the total profit
model.maximize(profit)
model.close()
# Parameterize the optimizer
optimizer.param.time_limit = time_limit
optimizer.solve()
#
# Write the solution
#
if output_file is not None:
with open(output_file, 'w') as f:
component_to_poduct = []
component_to_pool_fraction = []
pool_to_product = []
# Solution flows from the components to the products
for c in range(data.nbComponents):
for p in range(data.nbProducts):
component_to_poduct.append(
{"component": data.componentNames[c],
"product": data.productNames[p],
"flow": flowComponentToProduct[c][p].value})
# Solution fraction of the inflow at pool o coming from the component c
for c in range(data.nbComponents):
for o in range(data.nbPools):
component_to_pool_fraction.append(
{"component": data.componentNames[c],
"pool": data.poolNames[o],
"flow": fractionComponentToPool[c][o].value})
# Solution flows from the pools to the products
for o in range(data.nbPools):
for p in range(data.nbProducts):
pool_to_product.append(
{"pool": data.poolNames[o],
"product": data.productNames[p],
"flow": flowPoolToProduct[o][p].value})
json.dump({"objective": profit.value, "solution":
{"component_to_pool_fraction": component_to_pool_fraction,
"component_to_product": component_to_poduct,
"pool_to_product": pool_to_product}}, f)
if __name__ == '__main__':
if len(sys.argv) < 2:
print("Usage: python pooling.py instance_file [output_file] [time_limit]")
sys.exit(1)
instance_file = sys.argv[1]
output_file = sys.argv[2] if len(sys.argv) >= 3 else None
time_limit = int(sys.argv[3]) if len(sys.argv) >= 4 else 20
main(instance_file, output_file, time_limit)
- Compilation / Execution (Windows)
-
cl /EHsc pooling.cpp -I%HX_HOME%\include -Ilibs\lib-cpp /link %HX_HOME%\bin\hexaly130.libpooling instances\haverly1.json
- Compilation / Execution (Linux)
-
g++ pooling.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -Ilibs/lib-cpp -o pooling./pooling instances\haverly1.json
#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <json.hpp>
#include <map>
#include <sstream>
#include <vector>
using namespace hexaly;
using namespace std;
using json = nlohmann::json;
class Pooling {
public:
// Hexaly Optimizer
HexalyOptimizer optimizer;
// Instance data
int nbComponents;
int nbProducts;
int nbPools;
int nbAttributes;
vector<double> componentPrices;
vector<double> componentSupplies;
vector<vector<double>> componentQualities;
vector<string> componentNames;
map<string, int> componentIdx;
vector<double> productCapacities;
vector<double> productPrices;
vector<double> demand;
vector<vector<double>> maxTolerance;
vector<vector<double>> minTolerance;
vector<string> productNames;
map<string, int> productIdx;
vector<double> poolCapacities;
vector<string> poolNames;
map<string, int> poolIdx;
vector<vector<double>> upperBoundComponentToProduct;
vector<vector<double>> costComponentToProduct;
vector<vector<double>> upperBoundFractionComponentToPool;
vector<vector<double>> costComponentToPool;
vector<vector<double>> upperBoundPoolToProduct;
vector<vector<double>> costPoolToProduct;
vector<vector<vector<HxExpression>>> flowComponentToProductByPool;
// Decision variables
vector<vector<HxExpression>> flowComponentToProduct;
vector<vector<HxExpression>> flowPoolToProduct;
vector<vector<HxExpression>> fractionComponentToPool;
// Objective value
HxExpression profit;
/* Read instance data */
void readInstance(const string& fileName) {
ifstream infile;
infile.exceptions(ifstream::failbit | ifstream::badbit);
infile.open(fileName.c_str());
json problem;
problem << infile;
infile.close();
nbComponents = problem["components"].size();
nbProducts = problem["products"].size();
nbPools = problem["pool_size"].size();
nbAttributes = problem["components"][0]["quality"].size();
// Components
componentPrices.resize(nbComponents);
componentSupplies.resize(nbComponents);
componentQualities.resize(nbComponents);
componentNames.resize(nbComponents);
for (int c = 0; c < nbComponents; ++c) {
componentPrices[c] = problem["components"][c]["price"];
componentSupplies[c] = problem["components"][c]["upper"];
componentNames[c] = problem["components"][c]["name"].get<string>();
componentIdx[componentNames[c]] = c;
int currentAttributeIdx = 0;
componentQualities[c].resize(nbAttributes);
json qualityComponent = problem["components"][c]["quality"];
for (json::iterator it = qualityComponent.begin(); it != qualityComponent.end(); ++it) {
componentQualities[c][currentAttributeIdx] = it.value();
currentAttributeIdx++;
}
}
// Final products (blendings)
productCapacities.resize(nbProducts);
productPrices.resize(nbProducts);
demand.resize(nbProducts);
maxTolerance.resize(nbProducts);
minTolerance.resize(nbProducts);
productNames.resize(nbProducts);
for (int p = 0; p < nbProducts; ++p) {
productCapacities[p] = problem["products"][p]["upper"];
productPrices[p] = problem["products"][p]["price"];
demand[p] = problem["products"][p]["lower"];
productNames[p] = problem["products"][p]["name"].get<string>();
productIdx[productNames[p]] = p;
maxTolerance[p].resize(nbAttributes);
int currentAttributeIdx = 0;
json qualityUpperProduct = problem["products"][p]["quality_upper"];
for (json::iterator it = qualityUpperProduct.begin(); it != qualityUpperProduct.end(); ++it) {
maxTolerance[p][currentAttributeIdx] = it.value();
currentAttributeIdx++;
}
minTolerance[p].resize(nbAttributes);
currentAttributeIdx = 0;
json qualityLowerProduct = problem["products"][p]["quality_lower"];
bool hasLower = (qualityLowerProduct != NULL);
for (json::iterator it = qualityLowerProduct.begin(); it != qualityLowerProduct.end(); ++it) {
if (hasLower)
minTolerance[p][currentAttributeIdx] = it.value();
else
minTolerance[p][currentAttributeIdx] = 0;
currentAttributeIdx++;
}
}
// Intermediate pools
poolNames.resize(nbPools);
poolCapacities.resize(nbPools);
int currentPoolIdx = 0;
for (json::iterator it = problem["pool_size"].begin(); it != problem["pool_size"].end(); ++it) {
poolCapacities[currentPoolIdx] = it.value();
poolNames[currentPoolIdx] = it.key();
poolIdx[poolNames[currentPoolIdx]] = currentPoolIdx;
currentPoolIdx++;
}
// Flow graph
// Edges from the components to the products
upperBoundComponentToProduct.resize(nbComponents);
costComponentToProduct.resize(nbComponents);
flowComponentToProduct.resize(nbComponents);
// Edges from the components to the pools
upperBoundFractionComponentToPool.resize(nbComponents);
costComponentToPool.resize(nbComponents);
fractionComponentToPool.resize(nbComponents);
// Edges from the pools to the products
upperBoundPoolToProduct.resize(nbPools);
costPoolToProduct.resize(nbPools);
flowPoolToProduct.resize(nbPools);
// Flow from the components to the products and passing by the pools
flowComponentToProductByPool.resize(nbComponents);
for (int c = 0; c < nbComponents; ++c) {
upperBoundComponentToProduct[c].resize(nbProducts);
costComponentToProduct[c].resize(nbProducts);
flowComponentToProduct[c].resize(nbProducts);
upperBoundFractionComponentToPool[c].resize(nbPools);
costComponentToPool[c].resize(nbPools);
fractionComponentToPool[c].resize(nbPools);
flowComponentToProductByPool[c].resize(nbPools);
for (int o = 0; o < nbPools; ++o)
flowComponentToProductByPool[c][o].resize(nbProducts);
}
for (int o = 0; o < nbPools; ++o) {
upperBoundPoolToProduct[o].resize(nbProducts);
costPoolToProduct[o].resize(nbProducts);
flowPoolToProduct[o].resize(nbProducts);
}
// Bound and cost on the edges
for (json& edge : problem["component_to_product_bound"]) {
upperBoundComponentToProduct[componentIdx[edge["component"]]][productIdx[edge["product"]]] = edge["bound"];
if (edge.size() > 3)
costComponentToProduct[componentIdx[edge["component"]]][productIdx[edge["product"]]] = edge["cost"];
}
for (json& edge : problem["component_to_pool_fraction"]) {
upperBoundFractionComponentToPool[componentIdx[edge["component"]]][poolIdx[edge["pool"]]] =
edge["fraction"];
if (edge.size() > 3)
costComponentToPool[componentIdx[edge["component"]]][poolIdx[edge["pool"]]] = edge["cost"];
}
for (json& edge : problem["pool_to_product_bound"]) {
upperBoundPoolToProduct[poolIdx[edge["pool"]]][productIdx[edge["product"]]] = edge["bound"];
if (edge.size() > 3)
costPoolToProduct[poolIdx[edge["pool"]]][productIdx[edge["product"]]] = edge["cost"];
}
}
// Declare the optimization model
void solve(int limit) {
HxModel model = optimizer.getModel();
// Decision variables
// Flow from the components to the products
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p)
flowComponentToProduct[c][p] = model.floatVar(0, upperBoundComponentToProduct[c][p]);
}
// Fraction of the total flow in pool o coming from the component c
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o)
fractionComponentToPool[c][o] = model.floatVar(0, upperBoundFractionComponentToPool[c][o]);
}
// Flow from the pools to the products
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p)
flowPoolToProduct[o][p] = model.floatVar(0, upperBoundPoolToProduct[o][p]);
}
// Flow from the components to the products and passing by the pools
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p)
flowComponentToProductByPool[c][o][p] = fractionComponentToPool[c][o] * flowPoolToProduct[o][p];
}
}
// Constraints
// Proportion
for (int o = 0; o < nbPools; ++o) {
HxExpression proportion = model.sum();
for (int c = 0; c < nbComponents; ++c)
proportion.addOperand(fractionComponentToPool[c][o]);
model.constraint(proportion == 1);
}
// Component supply
for (int c = 0; c < nbComponents; ++c) {
HxExpression flowToProducts = model.sum(flowComponentToProduct[c].begin(), flowComponentToProduct[c].end());
HxExpression flowToPools = model.sum();
for (int o = 0; o < nbPools; ++o) {
flowToPools.addOperand(
model.sum(flowComponentToProductByPool[c][o].begin(), flowComponentToProductByPool[c][o].end()));
}
HxExpression totalOutFlow = model.sum(flowToPools, flowToProducts);
model.constraint(totalOutFlow <= componentSupplies[c]);
}
// Pool capacity (bounds on edges)
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
HxExpression flowComponentToPool =
model.sum(flowComponentToProductByPool[c][o].begin(), flowComponentToProductByPool[c][o].end());
HxExpression edgeCapacity = model.prod(poolCapacities[o], fractionComponentToPool[c][o]);
model.constraint(flowComponentToPool <= edgeCapacity);
}
}
// Product capacity
for (int p = 0; p < nbProducts; ++p) {
HxExpression flowFromPools = model.sum();
for (int o = 0; o < nbPools; ++o)
flowFromPools.addOperand(flowPoolToProduct[o][p]);
HxExpression flowFromComponents = model.sum();
for (int c = 0; c < nbComponents; ++c)
flowFromComponents.addOperand(flowComponentToProduct[c][p]);
HxExpression totalInFlow = model.sum(flowFromComponents, flowFromPools);
model.constraint(totalInFlow <= productCapacities[p]);
model.constraint(totalInFlow >= demand[p]);
}
// Product tolerance
for (int p = 0; p < nbProducts; ++p) {
for (int k = 0; k < nbAttributes; ++k) {
// Attribute from the components
HxExpression attributeFromComponents = model.sum();
for (int c = 0; c < nbComponents; ++c)
attributeFromComponents.addOperand(componentQualities[c][k] * flowComponentToProduct[c][p]);
// Attribute from the pools
HxExpression attributeFromPools = model.sum();
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o)
attributeFromPools.addOperand(componentQualities[c][k] * flowComponentToProductByPool[c][o][p]);
}
// Total flow in the blending
HxExpression totalFlowIn = model.sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowIn.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowIn.addOperand(flowPoolToProduct[o][p]);
HxExpression totalAttributeIn = model.sum(attributeFromComponents, attributeFromPools);
model.constraint(totalAttributeIn >= minTolerance[p][k] * totalFlowIn);
model.constraint(totalAttributeIn <= maxTolerance[p][k] * totalFlowIn);
}
}
// Objective function
// Cost of the flows from the components directly to the products
HxExpression directFlowCost = model.sum();
// Cost of the flows from the components to the products passing by the pools
HxExpression undirectFlowCost = model.sum();
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p) {
directFlowCost.addOperand(costComponentToProduct[c][p] * flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o) {
undirectFlowCost.addOperand((costComponentToPool[c][o] + costPoolToProduct[o][p]) *
flowComponentToProductByPool[c][o][p]);
}
}
}
// Gain of selling the final products
HxExpression productsGain = model.sum();
for (int p = 0; p < nbProducts; ++p) {
HxExpression totalFlowInProduct = model.sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowInProduct.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInProduct.addOperand(flowPoolToProduct[o][p]);
productsGain.addOperand(productPrices[p] * totalFlowInProduct);
}
// Cost of buying the components
HxExpression componentsCost = model.sum();
for (int c = 0; c < nbComponents; ++c) {
HxExpression totalFlowInComponent = model.sum();
for (int p = 0; p < nbProducts; ++p) {
totalFlowInComponent.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInComponent.addOperand(fractionComponentToPool[c][o] * flowPoolToProduct[o][p]);
}
componentsCost.addOperand(componentPrices[c] * totalFlowInComponent);
}
HxExpression totalCost = model.sum(componentsCost, directFlowCost, undirectFlowCost);
profit = model.sub(productsGain, totalCost);
// Maximize the total profit
model.maximize(profit);
model.close();
// Parameterize the optimizer
optimizer.getParam().setTimeLimit(limit);
optimizer.solve();
}
/* Write the solution */
void writeSolution(const string& fileName) {
ofstream outfile;
outfile.exceptions(ofstream::failbit | ofstream::badbit);
outfile.open(fileName.c_str());
vector<json> component_to_product;
vector<json> component_to_pool_fraction;
vector<json> pool_to_product;
component_to_product.resize(nbComponents * nbProducts);
component_to_pool_fraction.resize(nbComponents * nbPools);
pool_to_product.resize(nbPools * nbProducts);
// Solution flows from the components to the products
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p) {
component_to_product[c * nbProducts + p] = {{"component", componentNames[c]},
{"product", productNames[p]},
{"flow", flowComponentToProduct[c][p].getDoubleValue()}};
}
}
// Solution fraction of the inflow at pool o coming from the component c
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
component_to_pool_fraction[c * nbPools + o] = {
{"component", componentNames[c]},
{"pool", poolNames[o]},
{"flow", fractionComponentToPool[c][o].getDoubleValue()}};
}
}
// Solution flows from the pools to the products
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p) {
pool_to_product[o * nbProducts + p] = {{"pool", poolNames[o]},
{"product", productNames[p]},
{"flow", flowPoolToProduct[o][p].getDoubleValue()}};
}
}
json solution = {{"objective", profit.getDoubleValue()},
{"solution",
{{"component_to_pool_fraction", component_to_pool_fraction},
{"component_to_product", component_to_product},
{"pool_to_product", pool_to_product}}}};
outfile << solution;
}
};
int main(int argc, char** argv) {
if (argc < 2) {
cerr << "Usage: pooling inputFile [outputFile] [timeLimit] " << endl;
return 1;
}
const char* instanceFile = argv[1];
const char* solFile = argc > 2 ? argv[2] : NULL;
const char* strTimeLimit = argc > 3 ? argv[3] : "20";
try {
Pooling model;
model.readInstance(instanceFile);
model.solve(atoi(strTimeLimit));
if (solFile != NULL)
model.writeSolution(solFile);
return 0;
} catch (const exception& e) {
cerr << "An error occurred: " << e.what() << endl;
return 1;
}
}
- Compilation / Execution (Windows)
-
copy %HX_HOME%\bin\Hexaly.NET.dll .copy libs\lib-cs\Newtonsoft.Json.dll .csc Pooling.cs /reference:Hexaly.NET.dllPooling instances\haverly1.json
using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;
using Newtonsoft.Json.Linq;
using Newtonsoft.Json;
public class Pooling : IDisposable
{
// Hexaly Optimizer
private readonly HexalyOptimizer optimizer;
// Instance data
private int nbComponents;
private int nbProducts;
private int nbPools;
private int nbAttributes;
// Component data
private double[] componentPrices;
private double[] componentSupplies;
private double[,] componentQualities;
private string[] componentNames;
IDictionary<string, int> componentIdx;
// Final product data
private double[] productCapacities;
private double[] productPrices;
private double[] demand;
private double[,] maxTolerance;
private double[,] minTolerance;
private string[] productNames;
IDictionary<string, int> productIdx;
// Pool data
private double[] poolCapacities;
private string[] poolNames;
IDictionary<string, int> poolIdx;
// Flow graph data
private double[,] upperBoundComponentToProduct;
private double[,] costComponentToProduct;
private double[,] upperBoundFractionComponentToPool;
private double[,] costComponentToPool;
private double[,] upperBoundPoolToProduct;
private double[,] costPoolToProduct;
private HxExpression[,,] flowComponentToProductByPool;
// Decision variables
private HxExpression[,] flowComponentToProduct;
private HxExpression[,] fractionComponentToPool;
private HxExpression[,] flowPoolToProduct;
// Objective value
private HxExpression profit;
public Pooling()
{
optimizer = new HexalyOptimizer();
}
public void Dispose()
{
if (optimizer != null)
optimizer.Dispose();
}
// Read instance data
private void ReadInstance(string fileName)
{
JObject problem = JObject.Parse(File.ReadAllText(fileName));
// Components
JArray components = (JArray)problem["components"];
JObject attributes = (JObject)components[0]["quality"];
nbAttributes = attributes.Count;
nbComponents = components.Count;
componentPrices = new double[nbComponents];
componentSupplies = new double[nbComponents];
componentNames = new string[nbComponents];
componentIdx = new Dictionary<string, int>();
componentQualities = new double[nbComponents, nbAttributes];
int currentComponentIdx = 0;
foreach (JObject component in components)
{
componentPrices[currentComponentIdx] = (double)component["price"];
componentSupplies[currentComponentIdx] = (double)component["upper"];
componentNames[currentComponentIdx] = (string)component["name"];
componentIdx.Add(componentNames[currentComponentIdx], currentComponentIdx);
int currentAttributeIdx = 0;
attributes = (JObject)component["quality"];
for (
IEnumerator<KeyValuePair<String, JToken>> attributeIterator =
attributes.GetEnumerator();
attributeIterator.MoveNext();
)
{
KeyValuePair<String, JToken> attribute = attributeIterator.Current;
componentQualities[currentComponentIdx, currentAttributeIdx] =
(double)attribute.Value;
currentAttributeIdx++;
}
currentComponentIdx++;
}
// Products
JArray products = (JArray)problem["products"];
nbProducts = products.Count;
productPrices = new double[nbProducts];
productCapacities = new double[nbProducts];
demand = new double[nbProducts];
productNames = new string[nbProducts];
productIdx = new Dictionary<string, int>();
minTolerance = new double[nbProducts, nbAttributes];
maxTolerance = new double[nbProducts, nbAttributes];
int currentProductIdx = 0;
foreach (JObject product in products)
{
productPrices[currentProductIdx] = (double)product["price"];
productCapacities[currentProductIdx] = (double)product["upper"];
demand[currentProductIdx] = (double)product["lower"];
productNames[currentProductIdx] = (String)product["name"];
productIdx.Add(productNames[currentProductIdx], currentProductIdx);
int currentAttributeIdx = 0;
if (product["quality_lower"].Type != JTokenType.Null)
{
JObject qualityLower = (JObject)product["quality_lower"];
for (
IEnumerator<KeyValuePair<String, JToken>> qualityIterator =
qualityLower.GetEnumerator();
qualityIterator.MoveNext();
)
{
KeyValuePair<String, JToken> attribute = qualityIterator.Current;
minTolerance[currentProductIdx, currentAttributeIdx] = (double)attribute.Value;
currentAttributeIdx++;
}
}
currentAttributeIdx = 0;
JObject qualityUpper = (JObject)product["quality_upper"];
for (
IEnumerator<KeyValuePair<String, JToken>> qualityIterator =
qualityUpper.GetEnumerator();
qualityIterator.MoveNext();
)
{
KeyValuePair<String, JToken> attribute = qualityIterator.Current;
maxTolerance[currentProductIdx, currentAttributeIdx] = (double)attribute.Value;
currentAttributeIdx++;
}
currentProductIdx++;
}
// Intermediate pools
JObject pools = (JObject)problem["pool_size"];
nbPools = pools.Count;
poolCapacities = new double[nbPools];
poolNames = new string[nbPools];
poolIdx = new Dictionary<string, int>();
int currentPoolIdx = 0;
for (
IEnumerator<KeyValuePair<String, JToken>> poolIterator = pools.GetEnumerator();
poolIterator.MoveNext();
)
{
KeyValuePair<String, JToken> pool = poolIterator.Current;
poolCapacities[currentPoolIdx] = (double)pool.Value;
poolNames[currentPoolIdx] = (string)pool.Key;
poolIdx.Add(poolNames[currentPoolIdx], currentPoolIdx);
currentPoolIdx++;
}
// Flow graph
// Edges from the components to the products
upperBoundComponentToProduct = new double[nbComponents, nbProducts];
costComponentToProduct = new double[nbComponents, nbProducts];
flowComponentToProduct = new HxExpression[nbComponents, nbProducts];
// Edges from the components to the pools
upperBoundFractionComponentToPool = new double[nbComponents, nbPools];
costComponentToPool = new double[nbComponents, nbPools];
fractionComponentToPool = new HxExpression[nbComponents, nbPools];
// Edges from the pools to the products
upperBoundPoolToProduct = new double[nbPools, nbProducts];
costPoolToProduct = new double[nbPools, nbProducts];
flowPoolToProduct = new HxExpression[nbPools, nbProducts];
// Flow from the components to the products and passing by the pools
flowComponentToProductByPool = new HxExpression[nbComponents, nbPools, nbProducts];
// Bound and cost on the edges
foreach (JObject edge in problem["component_to_product_bound"].Value<JArray>())
{
string componentName = (String)edge["component"];
string productName = (String)edge["product"];
int componentIndex = (int)componentIdx[componentName];
int productIndex = (int)productIdx[productName];
double bound = (double)edge["bound"];
upperBoundComponentToProduct[componentIndex, productIndex] = bound;
if (edge.ContainsKey("cost"))
{
double cost = (double)edge["cost"];
costComponentToProduct[componentIndex, productIndex] = cost;
}
else
{
costComponentToProduct[componentIndex, productIndex] = 0;
}
}
foreach (JObject edge in problem["component_to_pool_fraction"].Value<JArray>())
{
string componentName = (string)edge["component"];
string poolName = (string)edge["pool"];
int componentIndex = (int)componentIdx[componentName];
int poolIndex = (int)poolIdx[poolName];
double bound = (double)edge["fraction"];
upperBoundFractionComponentToPool[componentIndex, poolIndex] = bound;
if (edge.ContainsKey("cost"))
{
double cost = (double)edge["cost"];
costComponentToPool[componentIndex, poolIndex] = cost;
}
else
{
costComponentToPool[componentIndex, poolIndex] = 0;
}
}
foreach (JObject edge in problem["pool_to_product_bound"].Value<JArray>())
{
string poolName = (string)edge["pool"];
string productName = (string)edge["product"];
int poolIndex = (int)poolIdx[poolName];
int productIndex = (int)productIdx[productName];
double bound = (double)edge["bound"];
upperBoundPoolToProduct[poolIndex, productIndex] = bound;
if (edge.ContainsKey("cost"))
{
double cost = (double)edge["cost"];
costPoolToProduct[poolIndex, productIndex] = cost;
}
else
{
costPoolToProduct[poolIndex, productIndex] = 0;
}
}
}
private void Solve(int limit)
{
// Declare the optimization model
HxModel model = optimizer.GetModel();
// Decision variables
// Flow from the components to the products
for (int c = 0; c < nbComponents; ++c)
{
for (int p = 0; p < nbProducts; ++p)
flowComponentToProduct[c, p] = model.Float(0, upperBoundComponentToProduct[c, p]);
}
// Fraction of the total flow in pool o coming from the component c
for (int c = 0; c < nbComponents; ++c)
{
for (int o = 0; o < nbPools; ++o)
fractionComponentToPool[c, o] = model.Float(
0,
upperBoundFractionComponentToPool[c, o]
);
}
// Flow from the pools to the products
for (int o = 0; o < nbPools; ++o)
{
for (int p = 0; p < nbProducts; ++p)
flowPoolToProduct[o, p] = model.Float(0, upperBoundPoolToProduct[o, p]);
}
// Flow from the components to the products and passing by the pools
for (int c = 0; c < nbComponents; ++c)
{
for (int o = 0; o < nbPools; ++o)
{
for (int p = 0; p < nbProducts; ++p)
flowComponentToProductByPool[c, o, p] = model.Prod(
fractionComponentToPool[c, o],
flowPoolToProduct[o, p]
);
}
}
// Constraints
// Proportion
for (int o = 0; o < nbPools; ++o)
{
HxExpression proportion = model.Sum();
for (int c = 0; c < nbComponents; ++c)
proportion.AddOperand(fractionComponentToPool[c, o]);
model.Constraint(model.Eq(proportion, 1));
}
// Component supply
for (int c = 0; c < nbComponents; ++c)
{
HxExpression flowToProducts = model.Sum();
HxExpression flowToPools = model.Sum();
for (int p = 0; p < nbProducts; ++p)
{
flowToProducts.AddOperand(flowComponentToProduct[c, p]);
for (int o = 0; o < nbPools; ++o)
flowToPools.AddOperand(flowComponentToProductByPool[c, o, p]);
}
HxExpression totalOutFlow = model.Sum(flowToPools, flowToProducts);
model.Constraint(model.Leq(totalOutFlow, componentSupplies[c]));
}
// Pool capacity (bounds on edges)
for (int c = 0; c < nbComponents; ++c)
{
for (int o = 0; o < nbPools; ++o)
{
HxExpression flowComponentToPool = model.Sum();
for (int p = 0; p < nbProducts; ++p)
flowComponentToPool.AddOperand(flowComponentToProductByPool[c, o, p]);
HxExpression edgeCapacity = model.Prod(
poolCapacities[o],
fractionComponentToPool[c, o]
);
model.Constraint(model.Leq(flowComponentToPool, edgeCapacity));
}
}
// Product capacity
for (int p = 0; p < nbProducts; ++p)
{
HxExpression flowFromPools = model.Sum();
for (int o = 0; o < nbPools; ++o)
flowFromPools.AddOperand(flowPoolToProduct[o, p]);
HxExpression flowFromComponents = model.Sum();
for (int c = 0; c < nbComponents; ++c)
flowFromComponents.AddOperand(flowComponentToProduct[c, p]);
HxExpression totalInFlow = model.Sum(flowFromComponents, flowFromPools);
model.Constraint(model.Leq(totalInFlow, productCapacities[p]));
model.Constraint(model.Geq(totalInFlow, demand[p]));
}
// Product tolerance
for (int p = 0; p < nbProducts; ++p)
{
for (int k = 0; k < nbAttributes; ++k)
{
// Attribute from the components
HxExpression attributeFromComponents = model.Sum();
for (int c = 0; c < nbComponents; ++c)
attributeFromComponents.AddOperand(
model.Prod(componentQualities[c, k], flowComponentToProduct[c, p])
);
// Attribute from the pools
HxExpression attributeFromPools = model.Sum();
for (int c = 0; c < nbComponents; ++c)
{
for (int o = 0; o < nbPools; ++o)
attributeFromPools.AddOperand(
model.Prod(
componentQualities[c, k],
flowComponentToProductByPool[c, o, p]
)
);
}
// Total flow in the blending
HxExpression totalFlowIn = model.Sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowIn.AddOperand(flowComponentToProduct[c, p]);
for (int o = 0; o < nbPools; ++o)
totalFlowIn.AddOperand(flowPoolToProduct[o, p]);
HxExpression totalAttributeIn = model.Sum(
attributeFromComponents,
attributeFromPools
);
model.Constraint(
model.Geq(totalAttributeIn, model.Prod(minTolerance[p, k], totalFlowIn))
);
model.Constraint(
model.Leq(totalAttributeIn, model.Prod(maxTolerance[p, k], totalFlowIn))
);
}
}
// Objective function
// Cost of the flows from the components directly to the products
HxExpression directFlowCost = model.Sum();
// Cost of the flows from the components to the products passing by the pools
HxExpression undirectFlowCost = model.Sum();
for (int c = 0; c < nbComponents; ++c)
{
for (int p = 0; p < nbProducts; ++p)
{
directFlowCost.AddOperand(
model.Prod(costComponentToProduct[c, p], flowComponentToProduct[c, p])
);
for (int o = 0; o < nbPools; ++o)
{
undirectFlowCost.AddOperand(
model.Prod(
costComponentToPool[c, o] + costPoolToProduct[o, p],
flowComponentToProductByPool[c, o, p]
)
);
}
}
}
// Gain of selling the final products
HxExpression productsGain = model.Sum();
for (int p = 0; p < nbProducts; ++p)
{
HxExpression totalFlowInProduct = model.Sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowInProduct.AddOperand(flowComponentToProduct[c, p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInProduct.AddOperand(flowPoolToProduct[o, p]);
productsGain.AddOperand(model.Prod(productPrices[p], totalFlowInProduct));
}
// Cost of buying the components
HxExpression componentsCost = model.Sum();
for (int c = 0; c < nbComponents; ++c)
{
HxExpression totalFlowInComponent = model.Sum();
for (int p = 0; p < nbProducts; ++p)
{
totalFlowInComponent.AddOperand(flowComponentToProduct[c, p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInComponent.AddOperand(
model.Prod(fractionComponentToPool[c, o], flowPoolToProduct[o, p])
);
}
componentsCost.AddOperand(model.Prod(componentPrices[c], totalFlowInComponent));
}
HxExpression totalCost = model.Sum(componentsCost, directFlowCost, undirectFlowCost);
profit = model.Sub(productsGain, totalCost);
// Maximize the total profit
model.Maximize(profit);
model.Close();
// Parameterize the optimizer
optimizer.GetParam().SetTimeLimit(limit);
optimizer.Solve();
}
/* Write the solution in a file */
private void WriteSolution(string fileName)
{
JArray component_to_product = new JArray();
JArray component_to_pool_fraction = new JArray();
JArray pool_to_product = new JArray();
// Solution flows from the components to the products
for (int c = 0; c < nbComponents; ++c)
{
for (int p = 0; p < nbProducts; ++p)
{
JObject solutionComponentToProduct = new JObject();
solutionComponentToProduct.Add("component", componentNames[c]);
solutionComponentToProduct.Add("product", productNames[p]);
solutionComponentToProduct.Add(
"flow",
flowComponentToProduct[c, p].GetDoubleValue()
);
component_to_product.Add(solutionComponentToProduct);
}
}
// Solution fraction of the inflow at pool o coming from the component c
for (int c = 0; c < nbComponents; ++c)
{
for (int o = 0; o < nbPools; ++o)
{
JObject solutionComponentToPoolFraction = new JObject();
solutionComponentToPoolFraction.Add("component", componentNames[c]);
solutionComponentToPoolFraction.Add("pool", poolNames[o]);
solutionComponentToPoolFraction.Add(
"flow",
fractionComponentToPool[c, o].GetDoubleValue()
);
component_to_pool_fraction.Add(solutionComponentToPoolFraction);
}
}
// Solution flows from the pools to the products
for (int o = 0; o < nbPools; ++o)
{
for (int p = 0; p < nbProducts; ++p)
{
JObject solutionPoolToProduct = new JObject();
solutionPoolToProduct = new JObject();
solutionPoolToProduct.Add("pool", poolNames[o]);
solutionPoolToProduct.Add("product", productNames[p]);
solutionPoolToProduct.Add("flow", flowPoolToProduct[o, p].GetDoubleValue());
pool_to_product.Add(solutionPoolToProduct);
}
}
JObject solution = new JObject();
solution.Add("component_to_product", component_to_product);
solution.Add("component_to_pool_fraction", component_to_pool_fraction);
solution.Add("pool_to_product", pool_to_product);
JObject output = new JObject();
output.Add("objective", profit.GetDoubleValue());
output.Add("solution", solution);
try
{
string outputString = JsonConvert.SerializeObject(output);
File.WriteAllText(fileName, outputString);
}
catch (Exception)
{
Console.WriteLine("Unable to write the solution");
}
}
public static void Main(string[] args)
{
if (args.Length < 1)
{
Console.WriteLine("Usage: Pooling inputFile [solFile] [timeLimit]");
Environment.Exit(1);
}
string instanceFile = args[0];
string outputFile = args.Length > 1 ? args[1] : null;
string strTimeLimit = args.Length > 2 ? args[2] : "20";
using (Pooling model = new Pooling())
{
model.ReadInstance(instanceFile);
model.Solve(int.Parse(strTimeLimit));
if (outputFile != null)
model.WriteSolution(outputFile);
}
}
}
- Compilation / Execution (Windows)
-
javac Pooling.java -cp %HX_HOME%\bin\hexaly.jar;libs\lib-java\gson-2.8.8.jarjava -cp %HX_HOME%\bin\hexaly.jar;libs\lib-java\gson-2.8.8.jar;. Pooling instances\haverly1.json
- Compilation / Execution (Linux)
-
javac Pooling.java -cp /opt/hexaly_13_0/bin/hexaly.jar;libs\lib-java\gson-2.8.8.jarjava -cp /opt/hexaly_13_0/bin/hexaly.jar:libs/lib-java/gson-2.8.8.jar:. Pooling instances/haverly1.json
import com.hexaly.optimizer.*;
import java.io.*;
import java.util.Map;
import com.google.gson.JsonObject;
import com.google.gson.JsonArray;
import com.google.gson.JsonParser;
import com.google.gson.stream.JsonReader;
public class Pooling {
// Hexaly Optimizer
private final HexalyOptimizer optimizer;
// Instance data
private int nbComponents;
private int nbProducts;
private int nbPools;
private int nbAttributes;
// Component data
private double[] componentPrices;
private double[] componentSupplies;
private double[][] componentQualities;
private String[] componentNames;
Map<String, Integer> componentIdx;
// Final product data
private double[] productCapacities;
private double[] productPrices;
private double[] demand;
private double[][] maxTolerance;
private double[][] minTolerance;
private String[] productNames;
Map<String, Integer> productIdx;
// Pool data
private double[] poolCapacities;
private String[] poolNames;
Map<String, Integer> poolIdx;
// Flow graph data
private double[][] upperBoundComponentToProduct;
private double[][] costComponentToProduct;
private double[][] upperBoundFractionComponentToPool;
private double[][] costComponentToPool;
private double[][] upperBoundPoolToProduct;
private double[][] costPoolToProduct;
private HxExpression[][][] flowComponentToProductByPool;
// Decision variables
private HxExpression[][] flowComponentToProduct;
private HxExpression[][] fractionComponentToPool;
private HxExpression[][] flowPoolToProduct;
// Objective value
private HxExpression profit;
private Pooling(HexalyOptimizer optimizer) {
this.optimizer = optimizer;
}
/* Read instance data */
private void readInstance(String fileName) throws IOException {
try {
JsonReader reader = new JsonReader(new InputStreamReader(new FileInputStream(fileName)));
JsonObject problem = JsonParser.parseReader(reader).getAsJsonObject();
// Components
JsonArray components = (JsonArray) problem.get("components");
nbComponents = components.size();
componentPrices = new double[nbComponents];
componentSupplies = new double[nbComponents];
componentNames = new String[nbComponents];
componentIdx = new java.util.HashMap<String, Integer>();
componentQualities = new double[nbComponents][];
int currentComponentIdx = 0;
for (Object objectComponent : components) {
JsonObject component = (JsonObject) objectComponent;
componentPrices[currentComponentIdx] = component.get("price").getAsDouble();
componentSupplies[currentComponentIdx] = component.get("upper").getAsDouble();
componentNames[currentComponentIdx] = component.get("name").getAsString();
componentIdx.put(componentNames[currentComponentIdx], currentComponentIdx);
int currentAttributeIdx = 0;
JsonObject attributes = (JsonObject) component.get("quality");
nbAttributes = attributes.size();
componentQualities[currentComponentIdx] = new double[nbAttributes];
for (Object objectKey : attributes.keySet()) {
String key = (String) objectKey;
componentQualities[currentComponentIdx][currentAttributeIdx] = attributes.get(key).getAsDouble();
currentAttributeIdx++;
}
currentComponentIdx++;
}
// Final products (blendings)
JsonArray products = (JsonArray) problem.get("products");
nbProducts = products.size();
productPrices = new double[nbProducts];
productCapacities = new double[nbProducts];
demand = new double[nbProducts];
productNames = new String[nbProducts];
productIdx = new java.util.HashMap<String, Integer>();
minTolerance = new double[nbProducts][nbAttributes];
maxTolerance = new double[nbProducts][nbAttributes];
int currentProductIdx = 0;
for (Object objectProduct : products) {
JsonObject product = (JsonObject) objectProduct;
productPrices[currentProductIdx] = product.get("price").getAsDouble();
productCapacities[currentProductIdx] = product.get("upper").getAsDouble();
demand[currentProductIdx] = product.get("lower").getAsDouble();
productNames[currentProductIdx] = product.get("name").getAsString();
productIdx.put(productNames[currentProductIdx], currentProductIdx);
int currentAttributeIdx = 0;
if (!product.get("quality_lower").isJsonNull()) {
JsonObject qualityLower = (JsonObject) product.get("quality_lower");
for (Object objectKey : qualityLower.keySet()) {
String key = (String) objectKey;
minTolerance[currentProductIdx][currentAttributeIdx] = qualityLower.get(key).getAsDouble();
currentAttributeIdx++;
}
}
currentAttributeIdx = 0;
JsonObject qualityUpper = (JsonObject) product.get("quality_upper");
for (Object objectKey : qualityUpper.keySet()) {
String key = (String) objectKey;
maxTolerance[currentProductIdx][currentAttributeIdx] = qualityUpper.get(key).getAsDouble();
currentAttributeIdx++;
}
currentProductIdx++;
}
// Intermediate pools
JsonObject pools = (JsonObject) problem.get("pool_size");
nbPools = pools.size();
poolCapacities = new double[nbPools];
poolNames = new String[nbPools];
poolIdx = new java.util.HashMap<String, Integer>();
int currentPoolIdx = 0;
for (Object objectKey : pools.keySet()) {
String key = (String) objectKey;
poolNames[currentPoolIdx] = key;
poolCapacities[currentPoolIdx] = pools.get(key).getAsDouble();
poolIdx.put(key, currentPoolIdx);
currentPoolIdx++;
}
// Flow graph
// Edges from the components to the products
upperBoundComponentToProduct = new double[nbComponents][nbProducts];
costComponentToProduct = new double[nbComponents][nbProducts];
flowComponentToProduct = new HxExpression[nbComponents][nbProducts];
// Edges from the components to the pools
upperBoundFractionComponentToPool = new double[nbComponents][nbPools];
costComponentToPool = new double[nbComponents][nbPools];
fractionComponentToPool = new HxExpression[nbComponents][nbPools];
// Edges from the pools to the products
upperBoundPoolToProduct = new double[nbPools][nbProducts];
costPoolToProduct = new double[nbPools][nbProducts];
flowPoolToProduct = new HxExpression[nbPools][nbProducts];
// Flow from the components to the products and passing by the pools
flowComponentToProductByPool = new HxExpression[nbComponents][nbPools][nbProducts];
// Bound and cost on the edges
for (Object objectEdge : (JsonArray) problem.get("component_to_product_bound")) {
JsonObject edge = (JsonObject) objectEdge;
String componentName = edge.get("component").getAsString();
String productName = edge.get("product").getAsString();
int componentIndex = (int) componentIdx.get(componentName);
int productIndex = (int) productIdx.get(productName);
double bound = edge.get("bound").getAsDouble();
upperBoundComponentToProduct[componentIndex][productIndex] = bound;
if (edge.has("cost")) {
double cost = edge.get("cost").getAsDouble();
costComponentToProduct[componentIndex][productIndex] = cost;
} else {
costComponentToProduct[componentIndex][productIndex] = 0;
}
}
for (Object objectEdge : (JsonArray) problem.get("component_to_pool_fraction")) {
JsonObject edge = (JsonObject) objectEdge;
String componentName = edge.get("component").getAsString();
String poolName = edge.get("pool").getAsString();
int componentIndex = (int) componentIdx.get(componentName);
int poolIndex = (int) poolIdx.get(poolName);
double bound = edge.get("fraction").getAsDouble();
upperBoundFractionComponentToPool[componentIndex][poolIndex] = bound;
if (edge.has("cost")) {
double cost = edge.get("cost").getAsDouble();
costComponentToPool[componentIndex][poolIndex] = cost;
} else {
costComponentToPool[componentIndex][poolIndex] = 0;
}
}
for (Object objectEdge : (JsonArray) problem.get("pool_to_product_bound")) {
JsonObject edge = (JsonObject) objectEdge;
String poolName = edge.get("pool").getAsString();
String productName = edge.get("product").getAsString();
int poolIndex = (int) poolIdx.get(poolName);
int productIndex = (int) productIdx.get(productName);
double bound = edge.get("bound").getAsDouble();
upperBoundPoolToProduct[poolIndex][productIndex] = bound;
if (edge.has("cost")) {
double cost = edge.get("cost").getAsDouble();
costPoolToProduct[poolIndex][productIndex] = cost;
} else {
costPoolToProduct[poolIndex][productIndex] = 0;
}
}
} catch (Exception ex) {
ex.printStackTrace();
}
}
private void solve(int limit) {
// Declare the optimization model
HxModel model = optimizer.getModel();
// Decision variables
// Flow from the components to the products
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p)
flowComponentToProduct[c][p] = model.floatVar(0, upperBoundComponentToProduct[c][p]);
}
// Fraction of the total flow in pool o coming from the component c
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o)
fractionComponentToPool[c][o] = model.floatVar(0, upperBoundFractionComponentToPool[c][o]);
}
// Flow from the pools to the products
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p)
flowPoolToProduct[o][p] = model.floatVar(0, upperBoundPoolToProduct[o][p]);
}
// Flow from the components to the products and passing by the pools
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p)
flowComponentToProductByPool[c][o][p] = model.prod(fractionComponentToPool[c][o],
flowPoolToProduct[o][p]);
}
}
// Constraints
// Proportion
for (int o = 0; o < nbPools; ++o) {
HxExpression proportion = model.sum();
for (int c = 0; c < nbComponents; ++c)
proportion.addOperand(fractionComponentToPool[c][o]);
model.constraint(model.eq(proportion, 1));
}
// Component supply
for (int c = 0; c < nbComponents; ++c) {
HxExpression flowToProducts = model.sum();
HxExpression flowToPools = model.sum();
for (int p = 0; p < nbProducts; ++p) {
flowToProducts.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
flowToPools.addOperand(flowComponentToProductByPool[c][o][p]);
}
HxExpression totalOutFlow = model.sum(flowToPools, flowToProducts);
model.constraint(model.leq(totalOutFlow, componentSupplies[c]));
}
// Pool capacity (bounds on edges)
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
HxExpression flowComponentToPool = model.sum();
for (int p = 0; p < nbProducts; ++p)
flowComponentToPool.addOperand(flowComponentToProductByPool[c][o][p]);
HxExpression edgeCapacity = model.prod(poolCapacities[o], fractionComponentToPool[c][o]);
model.constraint(model.leq(flowComponentToPool, edgeCapacity));
}
}
// Product capacity
for (int p = 0; p < nbProducts; ++p) {
HxExpression flowFromPools = model.sum();
for (int o = 0; o < nbPools; ++o)
flowFromPools.addOperand(flowPoolToProduct[o][p]);
HxExpression flowFromComponents = model.sum();
for (int c = 0; c < nbComponents; ++c)
flowFromComponents.addOperand(flowComponentToProduct[c][p]);
HxExpression totalInFlow = model.sum(flowFromComponents, flowFromPools);
model.constraint(model.leq(totalInFlow, productCapacities[p]));
model.constraint(model.geq(totalInFlow, demand[p]));
}
// Product tolerance
for (int p = 0; p < nbProducts; ++p) {
for (int currentAttributeIdx = 0; currentAttributeIdx < nbAttributes; ++currentAttributeIdx) {
// Attribute from the components
HxExpression attributeFromComponents = model.sum();
for (int c = 0; c < nbComponents; ++c)
attributeFromComponents.addOperand(
model.prod(componentQualities[c][currentAttributeIdx], flowComponentToProduct[c][p]));
// Attribute from the pools
HxExpression attributeFromPools = model.sum();
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o)
attributeFromPools.addOperand(model.prod(componentQualities[c][currentAttributeIdx],
flowComponentToProductByPool[c][o][p]));
}
// Total flow in the blending
HxExpression totalFlowIn = model.sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowIn.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowIn.addOperand(flowPoolToProduct[o][p]);
HxExpression totalAttributeIn = model.sum(attributeFromComponents, attributeFromPools);
model.constraint(
model.geq(totalAttributeIn, model.prod(minTolerance[p][currentAttributeIdx], totalFlowIn)));
model.constraint(
model.leq(totalAttributeIn, model.prod(maxTolerance[p][currentAttributeIdx], totalFlowIn)));
}
}
// Objective function
// Cost of the flows from the components directly to the products
HxExpression directFlowCost = model.sum();
// Cost of the flows from the components to the products passing by the pools
HxExpression undirectFlowCost = model.sum();
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p) {
directFlowCost.addOperand(model.prod(costComponentToProduct[c][p], flowComponentToProduct[c][p]));
for (int o = 0; o < nbPools; ++o) {
undirectFlowCost.addOperand(model.prod(costComponentToPool[c][o] + costPoolToProduct[o][p],
flowComponentToProductByPool[c][o][p]));
}
}
}
// Gain of selling the final products
HxExpression productsGain = model.sum();
for (int p = 0; p < nbProducts; ++p) {
HxExpression totalFlowInProduct = model.sum();
for (int c = 0; c < nbComponents; ++c)
totalFlowInProduct.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInProduct.addOperand(flowPoolToProduct[o][p]);
productsGain.addOperand(model.prod(productPrices[p], totalFlowInProduct));
}
// Cost of buying the components
HxExpression componentsCost = model.sum();
for (int c = 0; c < nbComponents; ++c) {
HxExpression totalFlowInComponent = model.sum();
for (int p = 0; p < nbProducts; ++p) {
totalFlowInComponent.addOperand(flowComponentToProduct[c][p]);
for (int o = 0; o < nbPools; ++o)
totalFlowInComponent.addOperand(model.prod(fractionComponentToPool[c][o], flowPoolToProduct[o][p]));
}
componentsCost.addOperand(model.prod(componentPrices[c], totalFlowInComponent));
}
HxExpression totalCost = model.sum(componentsCost, directFlowCost, undirectFlowCost);
profit = model.sub(productsGain, totalCost);
// Maximize the total profit
model.maximize(profit);
model.close();
// Parameterize the optimizer
optimizer.getParam().setTimeLimit(limit);
optimizer.solve();
}
/* Write the solution */
private void writeSolution(String fileName) throws IOException {
try {
JsonArray component_to_product = new JsonArray();
JsonArray component_to_pool_fraction = new JsonArray();
JsonArray pool_to_product = new JsonArray();
// Solution flows from the components to the products
for (int c = 0; c < nbComponents; ++c) {
for (int p = 0; p < nbProducts; ++p) {
JsonObject solutionComponentToProduct = new JsonObject();
solutionComponentToProduct.addProperty("component", componentNames[c]);
solutionComponentToProduct.addProperty("product", productNames[p]);
solutionComponentToProduct.addProperty("flow", flowComponentToProduct[c][p].getDoubleValue());
component_to_product.add(solutionComponentToProduct);
}
}
// Solution fraction of the inflow at pool o coming from the component c
for (int c = 0; c < nbComponents; ++c) {
for (int o = 0; o < nbPools; ++o) {
JsonObject solutionComponentToPoolFraction = new JsonObject();
solutionComponentToPoolFraction.addProperty("component", componentNames[c]);
solutionComponentToPoolFraction.addProperty("pool", poolNames[o]);
solutionComponentToPoolFraction.addProperty("flow", fractionComponentToPool[c][o].getDoubleValue());
component_to_pool_fraction.add(solutionComponentToPoolFraction);
}
}
// Solution flows from the pools to the products
for (int o = 0; o < nbPools; ++o) {
for (int p = 0; p < nbProducts; ++p) {
JsonObject solutionPoolToProduct = new JsonObject();
solutionPoolToProduct = new JsonObject();
solutionPoolToProduct.addProperty("pool", poolNames[o]);
solutionPoolToProduct.addProperty("product", productNames[p]);
solutionPoolToProduct.addProperty("flow", flowPoolToProduct[o][p].getDoubleValue());
pool_to_product.add(solutionPoolToProduct);
}
}
JsonObject solution = new JsonObject();
solution.add("component_to_product", component_to_product);
solution.add("component_to_pool_fraction", component_to_pool_fraction);
solution.add("pool_to_product", pool_to_product);
JsonObject output = new JsonObject();
output.addProperty("objective", profit.getDoubleValue());
output.add("solution", solution);
FileWriter file = new FileWriter(fileName);
file.write(output.toString().replaceAll("\\\\", ""));
file.close();
} catch (Exception ex) {
ex.printStackTrace();
}
}
public static void main(String[] args) {
if (args.length < 1) {
System.err.println("Usage: Pooling inaddFile [outaddFile] [timeLimit]");
System.exit(1);
}
String instanceFile = args[0];
String outaddFile = args.length > 1 ? args[1] : null;
String strTimeLimit = args.length > 2 ? args[2] : "20";
try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
Pooling model = new Pooling(optimizer);
model.readInstance(instanceFile);
model.solve(Integer.parseInt(strTimeLimit));
if (outaddFile != null) {
model.writeSolution(outaddFile);
}
} catch (Exception ex) {
System.err.println(ex);
ex.printStackTrace();
System.exit(1);
}
}
}