Pooling

Principles learned

  • Work with JSON files

  • Implement quadratic -bilinear- constraints

  • Use the operator sum

Problem

../_images/pooling.svg

In the pooling problem, flows of raw materials with specific attribute concentrations are blended to create final products whose attribute concentrations must be within tolerance intervals (attributes are interesting primary constituents of the materials). The flows can either be sent directly from the source nodes to the target nodes or can be blended at intermediate pools and then sent to the target nodes. 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 profit of the operation. The gains are the sale prices of the final products and the costs are the buying prices of the raw materials.

In that sense, the pooling problem can be seen as a mix between the minimal cost flow problem and the blending problem.

Download the example


Data

The instances come from the repository GitHub and are in format JSON. They are composed of the following elements:

  • “components” array of objects characterizing the raw materials with:

    • “name” name of the source node

    • “lower” minimal outflow (unused as equals to 0)

    • “upper” maximal outflow (supply)

    • “price” price of one unit of the raw material

    • “quality” object containing the concentration of each attribute in the raw material

  • “products” array of object characterizing the final products with:

    • “name” name of the target node

    • “lower” minimal inflow (demand)

    • “upper” maximal inflow (capacity)

    • “price” sale price of one unit of the final product

    • “quality_lower” object containing the minimal tolerate concentration of each attribute in the final product

    • “quality_upper” object containing the maximal tolerate 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 edge between the source nodes and the target nodes with:

    • “component” name of the source node c

    • “product” name of the target node p

    • “bound” maximal 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 edge between the source node and the pool nodes with:

    • “component” name of the source node c

    • “pool” name of the pool node o

    • “fraction” maximal 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 edge between the pool nodes and the target nodes with:

    • “pool” name of the pool node o

    • “product” name of the target node p

    • “bound” maximal 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 are either conveyed in the feature “price” or in the “cost” on the edges. The “cost” representation is used for the randstd instances and the “price” representation for the others.

We use JSON libraries for each API model to read the instance data and write the output solution: C# (Newtonsoft.Json), Java (gson-2.8.8), python (json), C++ (nlohmann/json.hpp). For Hexaly Modeler, the JSON module is used.

Program

This Hexaly Optimizer model is based on the Q formulation of the pooling problem which uses both flow proportions and flow quantities. The parameters of the problem are a set of raw materials C, a set of products P, a set of pools O and a tripartite flow graph composed with the three previous sets.

Three arrays of LSExpression decision variables are declared. They represent:

  • The flow from the source node c to the target node p for all (c, p) in CxP,

  • The flow from the pool o to the target node p for all (o, p) in OxP,

  • The proportion of the inflow at pool o coming from the source node c for all (c, o) in CxO.

An intermediate three-dimensional array of LSExpression is declared to represent the flow coming from the source node c and going to the target node p through the pool o which equals the proportion of inflow at pool o coming from the source node c times the flow from the pool o to the product p. Note that the flow from source c to pool o can be easily computed by summing the previous array over P for fixed c and o. This quantity is constrained by the proportion of inflow at pool o coming from the source c times the capacity of the pool.

For each pool, the total proportion of inflow coming from the source nodes are computed with the sum operator over all the components.This quantity is constrained to be equal to one.

The total inflow at the target nodes must satisfy the demand while respecting the capacity of the product. The total outflow at the source nodes cannot exceed the supply of the raw materials.

Finally, for each final product and each attribute, the concentration must be within a tolerance interval. This means the quantity of attribute in a final product must be within the interval times the total inflow at the target node. The quantity of attribute incoming at the target node is computed by summing the quantity coming directly from the source nodes and the quantity coming from the pools. The quantity coming from the source nodes is computed by summing over C the flow from a source node to the target node times the attribute concentration in the raw material. The quantity coming from the pools is computed by summing over CxO the flow coming from a source node to the target node through a pool times the concentration of attribute in the raw material.

The objective is to maximize the total profit of the operation. If the prices of the products and of the raw materials are explicit, the profit equals the sum over the products of the price times the total inflow minus the sum over the raw materials of the price times the total outflow. If the prices are expressed through costs on the edges, the profit is the opposite of the total cost. The total cost is computed by summing the unit cost times the flow over all the edges of the graph.

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\python
python pooling.py instances\haverly1.json
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python 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.lib
pooling 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.dll;libs\lib-cs\Newtonsoft.Json.dll
Pooling 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.jar
java -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.jar
java -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);
        }
    }
}