.. _sugarProduction:

Sugar production
================

The problem description in this section is taken from Section 6.4 *Cane sugar production* of the book *Applications of optimization with Xpress-MP*. The harvest of cane sugar in Australia is highly mechanized. The sugar cane is immediately transported to a sugarhouse in wagons that run on a network of small rail tracks. The sugar content of a wagonload depends on the field it has been harvested from and on the maturity of the sugar cane. Once harvested, the sugar content decreases rapidly through fermentation and the wagonload will entirely lose its value after a certain time. At this moment, eleven wagons loaded with the same quantity have arrived at the sugarhouse. They have been examined to find out the hourly loss and the remaining life span (in hours) of every wagon, these data are summarized in the following table:

+---------------+----+----+----+----+----+----+----+----+----+----+----+
| Lot           | 1  | 2  | 3  | 4  | 5  | 6  | 7  | 8  | 9  | 10 | 11 |
+---------------+----+----+----+----+----+----+----+----+----+----+----+
| Loss (kg/h)   | 43 | 26 | 37 | 28 | 13 | 54 | 62 | 49 | 19 | 28 | 30 |
+---------------+----+----+----+----+----+----+----+----+----+----+----+
| Life span (h) | 8  | 8  | 2  | 8  | 4  | 8  | 8  | 8  | 8  | 8  | 8  |
+---------------+----+----+----+----+----+----+----+----+----+----+----+

Every lot may be processed by any of the three, fully equivalent production lines of the sugarhouse. The processing of a lot takes two hours. It must be finished at the latest at the end of the life span of the wagonload. The manager of the sugarhouse wishes to determine a production schedule for the currently available lots that minimizes the total loss of sugar.

Model formulation
-----------------
Let :math:`WAGONS = {1, ..., NW}` be the set of wagons, :math:`NL` the number of production lines and :math:`DUR` the duration of the production process for every lot. The hourly loss for every wagon :math:`w` is given by :math:`LOSS_w` and its life span by :math:`LIFE_w`. We observe that, in an optimal solution, the production lines need to work without any break otherwise we could reduce the loss in sugar by advancing the start of the lot that follows the break. This means that the completion time of every lot is of the form :math:`s \times DUR`, with :math:`s > 0` and is an integer. The maximum value of :math:`s` is the number of time slots (of length :math:`DUR`) that the sugarhouse will work, namely :math:`NS = \text{ceil}(NW / NL)`, where *ceil* stands for *rounded to the next largest integer*. If :math:`NW / NL` is an integer, every line will process exactly :math:`NS` lots. Otherwise, some lines will process :math:`NS - 1` lots, but at least one line processes :math:`NS` lots. In all cases, the length of the optimal schedule is :math:`NS \times DUR` hours. We call :math:`SLOTS = \{1 ,... ,NS\}` the set of time slots. Every lot needs to be assigned to a time slot. We define variables :math:`process_w` for the time slot assigned to wagon :math:`w` and variables :math:`loss_w` for the loss incurred by this wagonload. Every time slot may take up to :math:`NL` lots because there are :math:`NL` parallel lines; therefore, we limit the number of *occurrences* of time slot values among the :math:`process_w` variables (this constraint relation is often called *cardinality constraint*) :

.. math::

    &\forall s \in SLOTS : |process_w = s|_{w \in WAGONS} \leq NL\\
    
The loss of sugar per wagonload :math:`w` and time slot :math:`s` is :math:`COST_{w,s} = s \times DUR \times LOSS_w`. Let variables :math:`loss_w` denote the loss incurred by wagon load :math:`w`:

.. math::

    &\forall w \in WAGONS : loss_w = COST_{w,proces_w}\\
 
The objective function (total loss of sugar) is then given as the sum of all losses:

.. math::

    &\text{minimize }\sum_{w \in WAGONS} loss_w\\

Implementation
-----------------

The following model is the implementation of this problem. It uses the function ceil to calculate the maximum number of time slots. The constraints on the processing variables are expressed by occurrence relations and the losses are obtained via element constraints. The branching strategy uses the variable selection criterion ``KLargestMin``, that is, choosing the variable with the largest lower bound.

.. tabs::
    
    .. code-tab:: c++
        
        // Number of wagon loads of sugar
        int NW = 11;
        // Number of production lines
        int NL = 3;
        // Time slots for production
        int NS = ceil(NW/((float)NL));
        // Loss in kg/hour
        int LOSS[] = { 43 ,26, 37, 28, 13, 54, 62, 49, 19, 28 ,30};
        // Remaining time per lot (in hours)
        int LIFE[] =  { 8 , 8,  2,  8 , 4 , 8 , 8,  8,  6,  8 , 8};    
        // Cost per wagon
        KIntArray COST;
        // Duration of the production
        int DUR = 2;
            
        // Loss per wagon
        KIntVarArray loss;  
            
        // Time slots for wagon loads
        KIntVarArray process;
            
        // Objective variable
        KIntVar * totalLoss;
            
        // Creation of the problem in this session
        KProblem problem(session,"A-4 Cane sugar production 1");

        int w;
        // Variables creation
        char name[80];
        for (w=0;w<NW;w++) {
            sprintf(name,"process(%i)",w);
            process += (* new KIntVar(problem,name,0,NS-1) ) ;        
        }
        for (w=0;w<NW;w++) {
            sprintf(name,"loss(%i)",w);
            loss += (* new KIntVar(problem,name,0,10000) ) ;
        }
            
        int s;
        // Wagon loads per time slot
        for (s=0; s < NS; s++) {
            KOccurTerm oc(s, process);        
            problem.post(new KOccurrence(oc, NL, false, true));            
        }
            
        // Limit on raw product life
        for (w = 0; w < NW; w++) {
            process[w].setSup(floor(LIFE[w]/((float)DUR))-1);
        }

        // initialization of COST array
        for (s = 0; s < NS; s++) {
            COST += 0;
        }
            
        KLinTerm totLossTerm;
        // Objective function: total loss  
        for (w = 0; w < NW; w++) {
            for (s=0;s<NS;s++) {
                COST[s] = (s+1)*DUR*LOSS[w];
            }
            KEltTerm kelt(COST,process[w]);
            problem.post(new KElement(kelt, loss[w],"element"));                
            totLossTerm = totLossTerm + loss[w];
        }    
        totalLoss = new KIntVar(problem,"totalLoss",0,10000);
        problem.post(totLossTerm == *totalLoss);
                
        // propagating problem
        if (problem.propagate()) {
            printf("Problem is infeasible\n");
            exit(1);
        }
            
        // search strategy customization
        KBranchingSchemeArray myBa;
        myBa += KAssignVar(KSmallestMax(),KMinToMax(),process);    
            
        // Solve the problem
        // creation of the solver
        KSolver solver(problem,myBa);
            
        // setting objective and sense of optimization
        problem.setSense(KProblem::Minimize);
        problem.setObjective(*totalLoss);
            
        int result = solver.optimize();
        // solution printing                
        KSolution * sol = &problem.getSolution();
        // print solution resume
        sol->printResume();
            
        // Solution printing
        printf("Total loss: %i\n", sol->getValue(*totalLoss));
        for (s = 0; s < NS; s++) {
            printf("Slot %i: ", s);
            for (w=0;w<NW;w++) {
                if (sol->getValue(process[w]) == s) {            
                    printf("wagon %i (%i) ",w,(s+1)*DUR*LOSS[w]);
                }
            }
            printf("\n");
        }
    
    .. code-tab:: py
    
        from kalis import *
        import math

        ### Data 
        nb_wagons = 11
        nb_production_lines = 3
        nb_time_slots = math.ceil(nb_wagons / float(nb_production_lines))

        # Loss in kg/hour
        wagons_loss_rate = [43, 26, 37, 28, 13, 54, 62, 49, 19, 28, 30]
        wagons_life_span = [8, 8, 2, 8, 4, 8, 8, 8, 6, 8, 8]

        # Duration of the production
        production_duration = 2

        ### Creation of the problem
        # Creation of the Kalis session
        session = KSession()

        # Creation of the optimization problem
        problem = KProblem(session, "A-4 Cane sugar production 1")

        ### Variables creation

        # Loss per wagon
        wagon_loss = KIntVarArray()
        for w in range(nb_wagons):
            wagon_loss += KIntVar(problem, "loss(%d)" % w, 0 , 10000)

        # Time slots for wagon loads
        wagon_process = KIntVarArray()
        for w in range(nb_wagons):
            wagon_process += KIntVar(problem, "process(%d)" % w, 0, nb_time_slots - 1)

        ### Set constraints
        # Set the cardinality constraint
        for s in range(nb_time_slots):
            oc = KOccurTerm(s, wagon_process)
            problem.post(KOccurrence(oc, nb_production_lines, False, True))

        # limit on raw product life:
        for w in range(nb_wagons):
            wagon_process[w].setSup(math.floor(wagons_life_span[w] / float(production_duration)) - 1)

        # Set the each wagon loss according to its processing order
        for w in range(nb_wagons):
            # Set cost per wagon
            cost = KIntArray()
            for t in range(nb_time_slots):
                res = cost.add( (t + 1) * production_duration * wagons_loss_rate[w] )
            # Add KElement constraint
            kelt = KEltTerm(cost, wagon_process[w])
            problem.post(KElement(kelt, wagon_loss[w], "element(%d)" % w))

        # Objective value
        wagons_loss_sum = 0
        for w in range(nb_wagons):
            wagons_loss_sum += wagon_loss[w]
        total_loss = KIntVar(problem, "totalLoss", 0, 10000)
        problem.post(total_loss == wagons_loss_sum)

        # First propagation to check inconsistency
        if problem.propagate():
            print("Problem is infeasible")
            sys.exit(1)

        ### Solve the problem
        # Search strategy customization
        my_branching_scheme = KBranchingSchemeArray()
        my_branching_scheme += KAssignVar(KSmallestMax(), KMinToMax(), wagon_process)

        # Set the solver
        solver = KSolver(problem, my_branching_scheme)

        # Setting objective and sense of optimization
        problem.setSense(KProblem.Minimize)
        problem.setObjective(total_loss)

        # Run optimization
        result = solver.optimize()
        if result:
            solution = problem.getSolution()
            solution.printResume()
            print("Total loss:", solution.getValue(total_loss))
            for t in range(nb_time_slots):
                print("Slot %d:" % t, end='')
                for w in range(nb_wagons):
                    if solution.getValue(wagon_process[w]) == t:
                        print(" wagon %d (%d)" % (w, (t + 1) * production_duration * wagons_loss_rate[w]), end='')
                print("")

    .. code-tab:: java
    
        // Number of wagon loads of sugar
        int NW = 11;
        // Number of production lines
        int NL = 3;
        // Time slots for production
        int NS = (int) Math.ceil(NW / (float) NL);
        // Loss in kg/hour
        int[] LOSS = {43, 26, 37, 28, 13, 54, 62, 49, 19, 28, 30};
        // Remaining time per lot (in hours)
        int[] LIFE = {8, 8, 2, 8, 4, 8, 8, 8, 6, 8, 8};
        // Cost per wagon
        KIntArray COST = new KIntArray();
        // Duration of the production
        int DUR = 2;

        // Loss per wagon
        KIntVarArray loss = new KIntVarArray();

        // Time slots for wagon loads
        KIntVarArray process = new KIntVarArray();

        // Objective variable
        KIntVar totalLoss = new KIntVar();

        // Creation of the session
        KSession session = new KSession();

        // Creation of the problem in this session
        KProblem problem = new KProblem(session, "A-4 Cane sugar production 1");

        int w;
        // Variables creation
        for (w = 0; w < NW; w++)
        {
            process.add(new KIntVar(problem, "process(" + w + ")", 0, NS - 1));
        }
        for (w = 0; w < NW; w++)
        {
            loss.add(new KIntVar(problem, "loss(" + w + ")", 0, 10000));
        }

        int s;
        // Wagon loads per time slot
        for (s = 0; s < NS; s++)
        {
            KOccurTerm oc = new KOccurTerm(s, process);
            problem.post(new KOccurrence(oc, NL, false, true));
        }

        // Limit on raw product life
        for (w = 0; w < NW; w++)
        {
            process.getElt(w).setSup(Math.floor(LIFE[w] / DUR) - 1);
        }

        KLinTerm totLossTerm = new KLinTerm();
        // Objective function: total loss
        for (w = 0; w < NW; w++)
        {
            for (s = 0; s < NS; s++)
            {
                COST.add((s + 1) * DUR * LOSS[w]);
            }
            KEltTerm kelt = new KEltTerm(COST,process.getElt(w));
            problem.post(new KElement(kelt, loss.getElt(w),"element"));
            totLossTerm.add(loss.getElt(w),1);
        }
        totalLoss = new KIntVar(problem,"totalLoss",0,10000);

        KNumVarArray intVarArrayToSet = totLossTerm.getLvars();
        KDoubleArray coeffsToSet = totLossTerm.getCoeffs();
        intVarArrayToSet.add(totalLoss);
        coeffsToSet.add(-1);
        problem.post(new KNumLinComb("",coeffsToSet,intVarArrayToSet,0,KNumLinComb.LinCombOperator.Equal));

        // propagating problem
        if (problem.propagate())
        {
            System.out.println("Problem is infeasible");
            exit(1);
        }

        // search strategy customization
        KBranchingSchemeArray myBa = new KBranchingSchemeArray();
        myBa.add(new KAssignVar(new KSmallestMax(), new KMinToMax(), process));

        // Solve the problem
        // creation of the solver
        KSolver solver = new KSolver(problem, myBa);

        // setting objective and sense of optimization
        problem.setSense(KProblem.Sense.Minimize);
        problem.setObjective(totalLoss);

        int result = solver.optimize();
        // solution printing
        KSolution sol = problem.getSolution();
        // print solution resume
        sol.printResume();

        // Solution printing
        System.out.println("Total loss: " + sol.getValue(totalLoss));
        for (s = 0; s < NS; s++)
        {
            System.out.print("Slot " + s + " :");
            for (w = 0; w < NW; w++)
            {
                if (sol.getValue(process.getElt(w)) == s)
                {
                    System.out.print("wagon " + w + " (" + (s + 1) * DUR * LOSS[w] + ")");
                }
            }
            System.out.println();
        }

An alternative formulation of the constraints on the processing variables is to replace the ``KOccurrence`` constraints by a single ``KGlobalCardinalityConstraint``, indicating for every time slot the minimum and maximum number (:math:`MINUSE_s = 0` and :math:`MAXUSEs = NL`) of production lines that may be used.

.. tabs::
    
    .. code-tab:: c++
        
        int * vals = new int[NS];
        int * low = new int[NS];
        int * upp = new int[NS];        
        for (s=0;s<NS;s++) {
            vals[s] = s;
            low[s] = 0;
            upp[s] = NL;
        }
        problem.post(new KGlobalCardinalityConstraint("Wagon_loads_per_time_slot\0", process, vals, low, upp,  NS));
        delete[] vals;
        delete[] low;
        delete[] upp;

    .. code-tab:: py
        
        vals = KIntArray()
        low = KIntArray()
        upp = KIntArray()
        for t in range(nb_time_slots):
            res = vals.add(t)
            res = low.add(0)
            res = upp.add(nb_production_lines)
        problem.post(KGlobalCardinalityConstraint("Wagon_loads_per_time_slot", wagon_process, vals, low, upp, nb_time_slots))

    .. code-tab:: java
    
        KIntArray vals = new KIntArray();
        KIntArray low = new KIntArray();
        KIntArray upp = new KIntArray();
        for (s = 0; s < NS; s++)
        {
            vals.add(s);
            low.add(0);
            upp.add(NL);
        }
        problem.post(new KGlobalCardinalityConstraint("Wagon_loads_per_time_slot", process, vals, low, upp, NS));

Yet another formulation of this problem is possible with Artelys Kalis, namely interpreting it as a cumulative scheduling problem where the wagon loads are represented by tasks of unit duration, scheduled on a discrete resource with a capacity corresponding to the number of production lines.

Results
-----------------

We obtain a total loss of 1620 kg of sugar. The corresponding schedule of lots is shown in the following table (there are several equivalent solutions) :

+----------------+----------------+----------------+----------------+
|    Slot 1      |    Slot 2      |    Slot 3      |    Slot 4      |
+----------------+----------------+----------------+----------------+
| lot 3 (74 kg)  | lot 1 (172 kg) | lot 4 (168 kg) | lot 2 (208 kg) |
+----------------+----------------+----------------+----------------+
| lot 6 (108 kg) | lot 5 (52 kg)  | lot 9 (114 kg) | lot 10 (224 kg)|
+----------------+----------------+----------------+----------------+
| lot 7 (124 kg) | lot 8 (196 kg) | lot 11 (180 kg)|                |
+----------------+----------------+----------------+----------------+