Pooling¶

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);
        }
    }
}