Pooling¶

Principles learned¶

Work with JSON files
Implement quadratic -bilinear- constraints

Use the operator sum

Problem¶

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 the LSP model, the JSON module is used.

Program¶

This localsolver 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.

LSP
Python
C++
C#
Java

Execution:: localsolver pooling.lsp inFileName=instances/haverly1.json [lsTimeLimit=] [solFileName=]

/********** pooling.lsp **********/

use json;

/* Read instance data */
function input() {
    local usage = "Usage: localsolver pooling.lsp "
        + "inFileName=inputFile [lsTimeLimit=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 
    componentPrices[c in 0..nbComponents-1] = problem["components"][c]["price"];
    componentSupplies[c in 0..nbComponents-1] = problem["components"][c]["upper"];
    componentQualities[c in 0..nbComponents-1] = problem["components"][c]["quality"].values();

    componentNames[c in 0..nbComponents-1] = problem["components"][c]["name"];
    for [c in 0..nbComponents-1]
        componentIdx[componentNames[c]] = c;

    // Final products (blendings)
    productPrices[p in 0..nbProducts-1] = problem["products"][p]["price"];
    productCapacities[p in 0..nbProducts-1] = problem["products"][p]["upper"];
    demand[p in 0..nbProducts-1] = problem["products"][p]["lower"];
    
    productNames[p in 0..nbProducts-1] = problem["products"][p]["name"];
    for [p in 0..nbProducts-1]
        productIdx[productNames[p]] = p;

    for [p in 0..nbProducts-1] {
        if (problem["products"][p]["quality_lower"] == NULL)
                minTolerance[p][k in 0..nbAttributes-1] = 0;
        else 
            minTolerance[p] = problem["products"][p]["quality_lower"].values();
    }
    maxTolerance[p in 0..nbProducts-1] = problem["products"][p]["quality_upper"].values();

    // Intermediate pools
    poolCapacities = problem["pool_size"].values();
    poolNames = problem["pool_size"].keys();
    for [o in 0..nbPools-1]
        poolIdx[poolNames[o]] = o;

    // Flow graph

    // Edges from the components to the products
    for [c in 0..nbComponents-1][p in 0..nbProducts-1] {
        upperBoundComponentToProduct[c][p] = 0;
        costComponentToProduct[c][p] = 0;
    }
    // Edges from the components to the pools
    for [c in 0..nbComponents-1][o in 0..nbPools-1] {
        upperBoundFractionComponentToPool[c][o] = 0;
        costComponentToPool[c][o] = 0;
    }
    // Edges from the pools to the products
    for [o in 0..nbPools-1][p in 0..nbProducts-1] {
        upperBoundPoolToProduct[o][p] = 0;
        costPoolToProduct[o][p] = 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
    for [c in 0..nbComponents-1][p in 0..nbProducts-1]
        flowComponentToProduct[c][p] <- float(0, upperBoundComponentToProduct[c][p]);

    // Fraction of the total flow in pool o coming from the component c
    for [c in 0..nbComponents-1][o in 0..nbPools-1]
        fractionComponentToPool[c][o] <- float(0, upperBoundFractionComponentToPool[c][o]);

    // Flow from the pools to the products
    for [o in 0..nbPools-1][p in 0..nbProducts-1]
        flowPoolToProduct[o][p] <- float(0, upperBoundPoolToProduct[o][p]);

    // Flow from the components to the products and passing by the pools
    for [c in 0..nbComponents-1][o in 0..nbPools-1][p in 0..nbProducts-1]
        flowComponentToProductByPool[c][o][p] <- fractionComponentToPool[c][o]
                * flowPoolToProduct[o][p];

    // Constraints

    // Proportion
    for [o in 0..nbPools-1] {
       proportion <- sum[c in 0..nbComponents-1](fractionComponentToPool[c][o]);
       constraint proportion == 1;
    }

    // Component supply
    for [c in 0..nbComponents-1] {
        flowToProducts <- sum[p in 0..nbProducts-1](flowComponentToProduct[c][p]);
        flowToPools <- sum[o in 0..nbPools-1][p in 0..nbProducts-1]
                (flowComponentToProductByPool[c][o][p]);
        totalOutFlow <- flowToPools + flowToProducts;
        constraint totalOutFlow <= componentSupplies[c];
    }

    // Pool capacity (bounds on edges)
    for [c in 0..nbComponents-1] {
        for [o in 0..nbPools-1] {
            flowComponentToPool <- sum[p in 0..nbProducts-1]
                    (flowComponentToProductByPool[c][o][p]);
            edgeCapacity <- poolCapacities[o] * fractionComponentToPool[c][o];
            constraint flowComponentToPool <= edgeCapacity;
        }
    }

    // Product capacity
    for [p in 0..nbProducts-1] {
        flowFromPools <- sum[o in 0..nbPools-1](flowPoolToProduct[o][p]);
        flowFromComponents <- sum[c in 0..nbComponents-1](flowComponentToProduct[c][p]);
        totalInFlow <- flowFromComponents + flowFromPools;
        constraint totalInFlow <= productCapacities[p];
        constraint totalInFlow >= demand[p];
    }

    // Product tolerance 
    for [p in 0..nbProducts-1] {
        for [k in 0..nbAttributes-1] {
            
            // Attribute from the components
            attributeFromComponents <- sum[c in 0..nbComponents-1]
                    (componentQualities[c][k] * flowComponentToProduct[c][p]);
            
            // Attribute from the pools
            attributeFromPools <- sum[c in 0..nbComponents-1](componentQualities[c][k] *
                    sum[o in 0..nbPools-1](flowComponentToProductByPool[c][o][p]));
            
            // Total flow in the blending
            totalFlowIn <- sum[o in 0..nbPools-1](flowPoolToProduct[o][p]) +
                    sum[c in 0..nbComponents-1](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-1][p in 0..nbProducts-1]
            (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-1][o in 0..nbPools-1][p in 0..nbProducts-1]
            ((costComponentToPool[c][o] + costPoolToProduct[o][p]) *
            flowComponentToProductByPool[c][o][p]);

    // Gain of selling the final products
    productsGain <- sum[p in 0..nbProducts-1](productPrices[p] *
            (sum[o in 0..nbPools-1](flowPoolToProduct[o][p]) +
            sum[c in 0..nbComponents-1](flowComponentToProduct[c][p])));
    
    // Cost of buying the components
    componentsCost <- sum[c in 0..nbComponents-1](componentPrices[c] *
            (sum[o in 0..nbPools-1][p in 0..nbProducts-1]
            (fractionComponentToPool[c][o] * flowPoolToProduct[o][p]) +
            sum[p in 0..nbProducts-1](flowComponentToProduct[c][p])));

    profit <- productsGain - componentsCost - (directFlowCost + undirectFlowCost);
    
    // Maximize the total profit
    maximize profit;
}

/* Parameterize the solver */
function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 20;
}

/* Write the solution */
function output() {

    if (solFileName == nil) return;

    // Solution flow from the components to the products
    for [c in 0..nbComponents-1] {
        for [p in 0..nbProducts-1] {
            component_to_poduct[c * nbProducts + p] = {"component" = componentNames[c],
                    "product" = productNames[p],
                    "flow" = flowComponentToProduct[c][p].value};
        }
    }
    // Solution fraction of the inflow at pool o coming from the component c
    for [c in 0..nbComponents-1] {
        for [o in 0..nbPools-1] {
            component_to_pool_fraction[c * nbPools + o] = {"component" = componentNames[c],
                    "pool" = poolNames[o],
                    "flow" = fractionComponentToPool[c][o].value};
        }
    }
    // Solution flows from the pools to the products
    for [o in 0..nbPools-1] {
        for [p in 0..nbProducts-1] {
            pool_to_product[o * nbProducts + p] = {"pool" = poolNames[o],
                    "product" = 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}}, solFileName);
}

Execution (Windows): set PYTHONPATH=%LS_HOME%\bin\python

python pooling.py instances\haverly1.json

Execution (Linux): export PYTHONPATH=/opt/localsolver_10_5/bin/python

python pooling.py instances/haverly1.json

########## pooling.py ##########

import localsolver
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 k 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 localsolver.LocalSolver() as ls:
        # Declare the optimization model
        model = ls.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 solver
        #
        if len(sys.argv) >= 4:
            ls.param.time_limit = int(sys.argv[3])
        else:
            ls.param.time_limit = 20

        ls.solve()

        # Write the solution
        if len(sys.argv) >= 3:
            with open(sys.argv[2], '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%LS_HOME%\include -Ilibs /link %LS_HOME%\bin\localsolver105.lib

pooling instances\haverly1.json

Compilation / Execution (Linux): g++ pooling.cpp -I/opt/localsolver_10_5/include -llocalsolver105 -lpthread -Ilibs -o pooling

./pooling instances/haverly1.json

//********* pooling.cpp *********

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <json.hpp>
#include <map>
#include "localsolver.h"

using namespace localsolver;
using namespace std;
using json = nlohmann::json;


class Pooling {
    // Solver 
    LocalSolver localsolver;

    // 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<LSExpression>>> flowComponentToProductByPool;

    // Decision variables
    vector<vector<LSExpression>> flowComponentToProduct;
    vector<vector<LSExpression>> flowPoolToProduct;
    vector<vector<LSExpression>> fractionComponentToPool;

    // Objective value
    LSExpression profit;

public:
    /* 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) {
        LSModel model = localsolver.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++) {
            LSExpression 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++) {
            LSExpression flowToProducts = model.sum(flowComponentToProduct[c].begin(),
                    flowComponentToProduct[c].end());
            LSExpression flowToPools = model.sum();
            for(int o = 0; o < nbPools; o++) {
                flowToPools.addOperand(model.sum(flowComponentToProductByPool[c][o].begin(),
                        flowComponentToProductByPool[c][o].end()));
            }
            LSExpression 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++) {
                LSExpression flowComponentToPool = model.sum(
                        flowComponentToProductByPool[c][o].begin(),
                        flowComponentToProductByPool[c][o].end());
                LSExpression edgeCapacity = model.prod(poolCapacities[o],
                        fractionComponentToPool[c][o]);
                model.constraint( flowComponentToPool <= edgeCapacity);
            }
        }

        // Product capacity
        for (int p = 0; p < nbProducts; p++) {
            LSExpression flowFromPools = model.sum();
            for (int o = 0; o < nbPools; o++) 
                flowFromPools.addOperand(flowPoolToProduct[o][p]);
            LSExpression flowFromComponents = model.sum();    
            for (int c = 0; c < nbComponents; c++)
                flowFromComponents.addOperand(flowComponentToProduct[c][p]);
            LSExpression 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
                LSExpression attributeFromComponents = model.sum();
                for (int c = 0; c < nbComponents; c++)
                    attributeFromComponents.addOperand(componentQualities[c][k] *
                            flowComponentToProduct[c][p]);

                // Attribute from the pools
                LSExpression 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
                LSExpression 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]);

                LSExpression 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
        LSExpression directFlowCost = model.sum();

        // Cost of the flows from the components to the products passing by the pools
        LSExpression 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
        LSExpression productsGain = model.sum();
        for (int p = 0; p < nbProducts; p++) {
            LSExpression 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
        LSExpression componentsCost = model.sum();
        for(int c = 0; c < nbComponents; c++) {
            LSExpression 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);
        }

        LSExpression totalCost = model.sum(componentsCost, directFlowCost, undirectFlowCost);
        profit = model.sub(productsGain, totalCost);

        // Maximize the total profit
        model.maximize(profit);
        model.close();

        // Parameterize the solver
        localsolver.getParam().setTimeLimit(limit);
        localsolver.getParam().setAdvancedParam("dumpFile", "dump_cpp.lsp");

        localsolver.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 %LS_HOME%\bin\localsolvernet.dll .

csc Pooling.cs /reference:localsolvernet.dll;libs\Newtonsoft.Json.dll

Pooling instances\haverly1.json

/********** Pooling.cs **********/

using System;
using System.IO;
using System.Collections.Generic;

using localsolver;
using Newtonsoft.Json.Linq;
using Newtonsoft.Json;

public class Pooling : IDisposable
{
    // Solver 
    private readonly LocalSolver localsolver;

    // 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 LSExpression[,,] flowComponentToProductByPool;

    // Decision variables
    private LSExpression[,] flowComponentToProduct;
    private LSExpression[,] fractionComponentToPool;
    private LSExpression[,] flowPoolToProduct;

    // Objective value
    private LSExpression profit;

    public Pooling()
    {
        localsolver = new LocalSolver();
    }

    public void Dispose()
    {
        if (localsolver != null)
            localsolver.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 LSExpression[nbComponents, nbProducts];

        // Edges from the components to the pools
        upperBoundFractionComponentToPool = new double[nbComponents, nbPools];
        costComponentToPool = new double[nbComponents, nbPools];
        fractionComponentToPool = new LSExpression[nbComponents, nbPools];

        // Edges from the pools to the products
        upperBoundPoolToProduct = new double[nbPools, nbProducts];
        costPoolToProduct = new double[nbPools, nbProducts];
        flowPoolToProduct = new LSExpression[nbPools, nbProducts];

        // Flow from the components to the products and passing by the pools
        flowComponentToProductByPool = new LSExpression[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
        LSModel model = localsolver.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++)
        {
            LSExpression 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++)
        {
            LSExpression flowToProducts = model.Sum();
            LSExpression 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]);
            }
            LSExpression 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++)
            {
                LSExpression flowComponentToPool = model.Sum();
                for (int p = 0; p < nbProducts; p++)
                    flowComponentToPool.AddOperand(flowComponentToProductByPool[c, o, p]);
                LSExpression edgeCapacity =  model.Prod(poolCapacities[o],
                        fractionComponentToPool[c, o]);
                model.Constraint(model.Leq(flowComponentToPool, edgeCapacity));
            }
        }

        // Product capacity
        for (int p = 0; p < nbProducts; p++)
        {
            LSExpression flowFromPools = model.Sum();
            for (int o = 0; o < nbPools; o++) 
                flowFromPools.AddOperand(flowPoolToProduct[o, p]);
            LSExpression flowFromComponents = model.Sum();
            for (int c = 0; c < nbComponents; c++)
                flowFromComponents.AddOperand(flowComponentToProduct[c, p]);
            LSExpression 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
                LSExpression attributeFromComponents = model.Sum();
                for (int c = 0; c < nbComponents; c++)
                    attributeFromComponents.AddOperand(model.Prod(
                            componentQualities[c, k], flowComponentToProduct[c, p]));

                // Attribute from the pools
                LSExpression 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
                LSExpression 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]);
                
                LSExpression 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
        LSExpression directFlowCost = model.Sum();

        // Cost of the flows from the components to the products passing by the pools
        LSExpression 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
        LSExpression productsGain = model.Sum();
        for (int p = 0; p < nbProducts; p++)
        {
            LSExpression 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
        LSExpression componentsCost = model.Sum();
        for(int c = 0; c < nbComponents; c++)
        {
            LSExpression 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));
        }

        LSExpression totalCost = model.Sum(componentsCost, directFlowCost, undirectFlowCost);
        profit = model.Sub(productsGain, totalCost);

        // Maximize the total profit
        model.Maximize(profit);
        model.Close();

        // Parameterize the solver
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.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 %LS_HOME%\bin\localsolver.jar;libs\gson-2.8.8.jar

java -cp %LS_HOME%\bin\localsolver.jar;libs\gson-2.8.8.jar;. Pooling instances\haverly1.json

Compilation / Execution (Linux): javac Pooling.java -cp %LS_HOME%/bin/localsolver.jar:libs/gson-2.8.8.jar

java -cp %LS_HOME%/bin/localsolver.jar:libs/gson-2.8.8.jar:. Pooling instances/haverly1.json

/********** Pooling.java **********/

import localsolver.*;
import java.io.*;
import java.util.Map;
import com.google.gson.JsonObject;
import com.google.gson.JsonArray;
import com.google.gson.JsonElement;
import com.google.gson.JsonParser;
import com.google.gson.stream.JsonReader;

public class Pooling {
    // Solver
    private final LocalSolver localsolver;

    // 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 LSExpression[][][] flowComponentToProductByPool;

    // Decision variables
    private LSExpression[][] flowComponentToProduct;
    private LSExpression[][] fractionComponentToPool;
    private LSExpression[][] flowPoolToProduct;

    // Objective value
    private LSExpression profit;

    private Pooling(LocalSolver localsolver) {
        this.localsolver = localsolver;
    }

    /* 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 LSExpression[nbComponents][nbProducts];
            
            // Edges from the components to the pools
            upperBoundFractionComponentToPool = new double[nbComponents][nbPools];
            costComponentToPool = new double[nbComponents][nbPools];
            fractionComponentToPool = new LSExpression[nbComponents][nbPools];
            
            // Edges from the pools to the products
            upperBoundPoolToProduct = new double[nbPools][nbProducts];
            costPoolToProduct = new double[nbPools][nbProducts];
            flowPoolToProduct = new LSExpression[nbPools][nbProducts];

            // Flow from the components to the products and passing by the pools
            flowComponentToProductByPool = new LSExpression[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();
        }
    }

    /* Declare the optimization model */
    private void solve(int limit) {
        LSModel model = localsolver.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++) {
            LSExpression 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++) {
            LSExpression flowToProducts = model.sum();
            LSExpression 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]);
            }
            LSExpression 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++) {
                LSExpression flowComponentToPool = model.sum();
                for (int p = 0; p < nbProducts; p++)
                    flowComponentToPool.addOperand(flowComponentToProductByPool[c][o][p]);
                LSExpression edgeCapacity = model.prod(poolCapacities[o],
                        fractionComponentToPool[c][o]);
                model.constraint(model.leq(flowComponentToPool, edgeCapacity));
            }
        }

        // Product capacity
        for (int p = 0; p < nbProducts; p++) {
            LSExpression flowFromPools = model.sum();
            for (int o = 0; o < nbPools; o ++) 
                flowFromPools.addOperand(flowPoolToProduct[o][p]);
            LSExpression flowFromComponents = model.sum();
            for (int c = 0; c < nbComponents; c++)
                flowFromComponents.addOperand(flowComponentToProduct[c][p]);
            LSExpression 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
                LSExpression attributeFromComponents = model.sum();
                for (int c = 0; c < nbComponents; c++)
                    attributeFromComponents.addOperand(model.prod(
                            componentQualities[c][currentAttributeIdx],
                            flowComponentToProduct[c][p]));

                // Attribute from the pools
                LSExpression 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
                LSExpression 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]);

                LSExpression 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
        LSExpression directFlowCost = model.sum();
        // Cost of the flows from the components to the products passing by the pools
        LSExpression 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
        LSExpression productsGain = model.sum();
        for (int p = 0; p < nbProducts; p++) {
            LSExpression 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
        LSExpression componentsCost = model.sum();
        for(int c = 0; c < nbComponents; c++) {
            LSExpression 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));
        }

        LSExpression totalCost = model.sum(componentsCost, directFlowCost, undirectFlowCost);
        profit = model.sub(productsGain, totalCost);
        
        // Maximize the total profit
        model.maximize(profit);
        model.close();

        // Parameterize the solver
        localsolver.getParam().setTimeLimit(limit);
        
        localsolver.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 (LocalSolver localsolver = new LocalSolver()) {
            Pooling model = new Pooling(localsolver);
            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);
        }
    }
}