- Project survey
This project concerns production planning at three paper mills in Sweden. Briefly, the raw materials, assortments, are transported from different forest areas, districts, to the inventory of the mills. The assortments are then used in the production process to produce different products. Every product requires a certain amount of the different assortments, which must be available when the production of that product begins. Afterwards, the products are transported to the inventory before transportation to the international and domestic customers.
The task is to investigate if it is possible to aggregate the production at the mills and hopefully produce better solutions than when the production at each mill is optimized separately.
The concept of production plans
The production at the different mills is governed by recipes. The recipes state how much of each assortment that is needed and how much of each product that is produced. It is assumed that one and only one recipe is used every day, but different mills can use different recipes. As an example, to run RECIPE1 at VARO requires 20.8 units of GRAN and 19.2 units of LOV. The output of this recipe is 30 units of the product S70.
The production planning is performed for a time period of 30 days. In this project, two types of costs associated with the production are used. These are the cost of using a recipe and the cost of switching between recipes.
Basic data
Since this project only concerns the production, how the assortments are provided is irrelevant and will not be further discussed. However, in order to understand the process, the underlying data is presented as follows.
Districts:EAST, SOUTH, WEST, CENTRAL NORTH
Assortments:TAL, GRAN, LOV, ASP
Mills:VARO, MONSTRAS, MORRUM
Products:S70, S80, S85, S90, SBZ
Recipes:RECIPE1, RECIPE2, RECIPE3, RECIPE4, RECIPE5
Domestic customers:MALMO, GOTEBORG, STOCKHOLM
Overseas customers:UK1, UK2, GERMAN1, FRANCE1, FRANCE2, ITALY1, SPAIN1
The process in its entirety, from providing the assortments to the delivery, is depicted in figure 1.
The problem and solution method
The problem is outlined as followed. From an initial set of recipes, one for each mill, used on day zero, the problem is to define which set of recipes to use every day. The recipes should be used in such ways that the cost of production is minimized and the demand of all customers is satisfied.
Since the demands at the different customers are known in advance, but the transportation times are not, it is necessary to create a total demand that the production must satisfy. This was done by shifting the domestic and oversea demands. If the shift is one and three days respectively, this means that the domestic demands must be satisfied one day before they actually occur and the oversea demands three days before they actually occur. This is to simulate a transportation time of one and three days respectively.
The solution method is a well-known approach by the name column generation, which is described in section II. Thus, the problem is divided into two problems that are easier to solve. One of these problems, the subproblem, takes care of the optimization part and the other one, the master problem, takes care of the combinatorial part.
Figure 1 The process
- Column generation
Some of the linear programming problem encountered in practice may be very large in terms of the number of variables and/or constraints. Decomposition principles is one of the mostly used solution methods applied to such problems to obtain the solution. According to this principle, the original problem is decomposed into small subproblems and then these subproblems are solved independently. Often, the information created in one of the problems is used to guide the other problems.
The procedure, when applicable, has the advantage of making it possible to solve large-scale problems which otherwise may be computationally very difficult to solve.
The idea behind column generation is to solve a reduced version of the real problem and then use information from this solution to add essential variables. This procedure is then repeated until no more essential variables can be added.
In this project, the reduced real problem, the master problem, is a LP-relaxed version of the real problem. Since this approach uses the dual prices of different constraints, a LP-relaxation of the real problem is used. The real problem has access to all possible production plans, but the number of possible plans is too large for any solver to handle. In the reduced problem, only a subset of the production plans is used. The master problem is solved with this set of production plans, and the dual information is saved.
With the dual information of the master problem, a new problem, the subproblem, is solve to decide which of the production plans not in the master problem that should be added to the master problem. If such a production plan is found, it is added to the master problem, which is then re-solved.
This procedure is often called column generation procedure. The objective function of master problem expresses the total cost of the overall problem and the objective function of subproblem expresses how good the generated column is. Figure 2 depict the decomposition principle and switching of information between the master problem and the subproblem.
Figure 2 Decomposition principle
The dual prices of the master problem are used in the objective function of the subproblem to guide the solution procedure of the subproblem. Often there is an attractive structure in the subproblem, which makes the subproblem easy to solve. This means that particular attention should be paid to the different subproblems that can arise and to the variety of methods utilized to solve them. These include dynamic programming, linear programming, nonlinear programming and network flow algorithms.
There are a variety of termination criteria based on the accuracy of the desired solution to the overall problem. Two such criteria are as follows.
- UBD-LBD termination criteria
An upper bound (UBD) and lower bound (LBD) to the relaxed master problem is calculated and the program is terminated when the gap between these bounds is less than a given constant.
- Optimal termination criteria
This is the simplest criterion that can be used to terminate the execution of the program. The program is terminated when the subproblem does not find any new production plans to add to the master problem. The bottleneck of this criterion is that it can be very time consuming and hence is seldom applied in practice.
In this project non of the above criteria are used, but this will be explained later.
- Decomposition
The master problem
In this case, the master problem is the linear relaxation of the original problem. The objective function in master is the overall cost of the whole system over the time period of 30 days. These costs include:
- The total cost of transporting assortments to the pulp mills.
- The total cost of storage of assortments at the districts.
- The total cost of production plans.
- The total cost of transportation of products from pulp mills to the harbor ports.
- The total cost of using ships schedules.
- The total cost of shortage of products for domestic customers.
- The total cost of oversupply of products for domestic customers.
- The total cost of shortage of products for overseas customers.
- The total cost of oversupply of products for overseas customers.
- The total cost of storage of products at the harbor ports.
Formulation of the master problem is outlined in the file by the name model.mod and the reader is referred to this file for further considerations.
The subproblem
The objective function in the subproblem is the total reduced costs for all pulp mills. The reduced cost of each pulp mill is defined as the sum of the three different costs listed below.
- Cost of using a recipe.
- Cost of switching between recipes.
- Cost of using a production plan.
The costs of using a recipe and switching between recipes are fixed costs only depending on the production plan used and not on the solution of the master problem. The cost of using a production plan is calculated by using the dual prices obtained from solving the master problem. These dual prices (as defined in the program) are:
- FCA_PRICES
The dual variables corresponding to the flow conservation constraint for assortments.
- FCP_PRICES
The dual variables corresponding to the flow conservation constraint for products.
- CON_PRICES
The dual variables corresponding to the convexity constraint. This constraint expresses that one and only one production plan should be used on every pulp mill.
These flow conservation constraints for assortments and products and the convexity constraint are the only constraints in which the variables prodplanusage, and hence the only constraints of our concern.
With these dual prices it is possible to define the subproblem.
Definition of parameters
param FCA_PRICES {MILLS, ASSORTS, PERIODS} default 0.0;
The dual variables corresponding to the flow conservation constraint for assortments. There are 360 dual variables of this type.
param FCP_PRICES {MILLS, PRODS, PERIODS} default 0.0;
The dual variables corresponding to the flow conservation constraint for products. There are 450 variables of this type.
param CON_PRICES {MILLS} default 0.0;
The dual variables corresponding to the convexity constraint. There are 3 dual variables of this type.
param TOTAL_DEMAND{PERIODS, PRODS};
This parameter contains the total demand of domestic and international customers. It is calculated as: TOTAL_DEMAND[t,p] :=
sum {d in DDEMS} DEMAND_DDEMS[d,p,t+shift_domestic]) +
sum {c in CUSTS} DEMAND_CUST[c,p,t+shift_international]);
Definition of variables
var Y {m in MILLS, r in RECIPES[m], t in PERIODS} binary;
var X {m in MILLS, r1 in RECIPES[m], r2 in RECIPES[m], t in PERIODS} >= 0;
var inventory {t in PERIODS, p in PRODS} >= 0;
The inventory of product p at the end of day t
Objective function
Minimize total_reduced_cost:
The total cost of using recipes.
sum{m in MILLS, r in RECIPES[m], t in PERIODS} C_RECIPE[m, r] * Y[m, r, t]
+
The total cost of switching between recipes.
sum{m in MILLS, r1 in RECIPES[m], r2 in RECIPES[m], t in PERIODS}
C_RECIPECHANGE[m, r1, r2] * X[m, r1, r2, t]
-
The total cost corresponding to the relaxed flow conservation constraint for assortments.
sum{m in MILLS, r in RECIPES[m], a in ASSORTS, t in PERIODS}
RECIPE_INPUT[m, r, a] * Y[m, r, t] * FCA_PRICES[m, a, t]
+
The total cost corresponding to the relaxed flow conservation constraint for products.
sum{m in MILLS, r in RECIPES[m], p in PRODS, t in PERIODS}
RECIPE_OUTPUT[m, r, p] * Y[m, r, t] * FCP_PRICES[m, p, t]
-
The total cost corresponding to the relaxed convexity constraint.
sum{m in MILLS} CON_PRICES[m];
The constraints
The inventory constraints.
Subject to INVENTORY_BALANCE {t in PERIODS, p in PRODS}:
inventory[t,p] = sum { m in MILLS, r in RECIPES[m]} Y[m,r,t] * RECIPE_OUTPUT[m,r,p] +
if t = 1 then sum {m in MILLS} INIT_STORE_MILL_PROD[m,p] else inventory[t-1,p]) - TOTAL_DEMAND[t,p];
Constraint that ensures one recipe is used each day on each mill.
Subject to ONE_RECIPE_MILL_AND_DAY {m in MILLS, t in PERIODS}:
sum{r in RECIPES[m]}Y[m, r, t]= 1;
This constraint connects the variables X and Y for the initial period.
Subject to INITIAL_CONNECTION {m in MILLS, r in RECIPES[m]}:
Y[m,r,1] <= X[m,INIT_RECIPE[m],r,1];
These constraints connect the variables X and Y.
Subject to CONNECTION
{t in PERIODS, m in MILLS, r1 in RECIPES[m], r2 in RECIPES[m]: t > first(PERIODS)}:
Y[m,r1,t-1] + Y[m,r2,t] - 1 <= X[m,r1,r2,t];
The solution of the subproblem is a combined production plan including one production plan for each mill. This combined production plan is then divided before the master problem is solved.
The subproblem is an integer problem in nature, but the problem does not fall in any of the known problem classes.
The subproblem was solved by CPLEX with a branch-and-bound method. The termination criteria used for the master problem were a specific number of production plans generated and for the subproblem both a UBD-LBD criteria and a specific number of found solutions.
- Results
In the first runs (see table 1), 15 production plans were generated in the column generation process with different shifts in the domestic (DDEM shift) and the international (IDEM shift) demand. The termination criteria used for the subproblem in all the runs except for the last, were either 2 feasible solutions or a gap of less than 1 %. When the 15 plans were generated, the master problem was solved 15 times, one for each of the generated plans where the specific plan was fixed to one, in order to evaluate it.
Settings / ResultsRun / Plans / DDEM shift / IDEM shift / Term crit / Opt Master LP / Opt Master IP / With plan nr
1 / 15 / 0 / 3 / 2, 1% / 868617 / 895198 / 13
2 / 15 / 0 / 4 / 2, 1% / 853078 / 898153 / 12
3 / 15 / 0 / 5 / 2, 1% / 857914 / 897529 / 14
4 / 15 / 0 / 6 / 2, 1% / 852186 / 893695 / 13
5 / 15 / 1 / 3 / 2, 1% / 853582 / 898125 / 9
6 / 15 / 1 / 4 / 2, 1% / 853413 / 896557 / 14
7 / 15 / 1 / 5 / 2, 1% / 864057 / 902541 / 15
8 / 15 / 1 / 6 / 2, 1% / 857567 / 900217 / 7
9 / 15 / 2 / 3 / 2, 1% / 865120 / 895472 / 15
10 / 15 / 2 / 4 / 2, 1% / 847925 / 904265 / 9
11 / 15 / 2 / 5 / 2, 1% / 861890 / 891606 / 15
12 / 15 / 2 / 6 / 2, 1% / 864852 / 892584 / 13
14 / 15 / 1 / 4 / 3, 5% / 853778 / 896959 / 12
Table 1: Runs where the production plans were fixed when solving the master problem.
As we observe in table 1, there’s a large difference between the LP solution (Opt Master LP) for the master and the solution where the production plans are fixed (Opt Master IP). We felt that this difference would decrease if the different mills were not forced to use a plan from the same subproblem generation. Therefore, in the second set of runs (see table 2), instead of fixing which production plans that should be used, we only redeclared the prodplanusage variables to binary and solved the master problem as an integer problem.
Settings / ResultsRun / Plans / DDEM shift / IDEM shift / Term crit / Opt Master LP / Opt Master IP / With plan nr
Monsteras / Morrum / Varo
15 / 15 / 0 / 3 / 2, 1% / 868617 / 872452 / 14 / 15 / 19
16 / 15 / 0 / 4 / 2, 1% / 855234 / 855722 / 17 / 11 / 20
17 / 15 / 0 / 5 / 2, 1% / 857914 / 861305 / 19 / 12 / 14
18 / 15 / 0 / 6 / 2, 1% / 846323 / 846783 / 16 / 16 / 12
19 / 15 / 1 / 3 / 2, 1% / 853582 / 856574 / 20 / 14 / 14
20 / 15 / 1 / 4 / 2, 1% / 853413 / 856676 / 19 / 16 / 20
21 / 15 / 1 / 5 / 2, 1% / 864057 / 865140 / 10 / 14 / 15
22 / 15 / 1 / 6 / 2, 1% / 857567 / 859812 / 12 / 18 / 13
23 / 15 / 2 / 3 / 2, 1% / 865120 / 866688 / 15 / 11 / 9
24 / 15 / 2 / 4 / 2, 1% / 847925 / 851195 / 11 / 13 / 19
25 / 15 / 2 / 5 / 2, 1% / 861890 / 865655 / 14 / 9 / 20
26 / 15 / 2 / 6 / 2, 1% / 864852 / 865910 / 10 / 14 / 15
Table 2: Runs where the master was solved as an IP.
In table 2, the last three columns tell which production plan each of the mills use in the optimal solution. As there were five (six in one case) plans to start which and we generated 15 more, it is possible to use plan number 20 as Varo does in run 16. The settings for the subproblem were the same as in table 1, but the results are quite more pleasing. The LP solutions are of course still better than the IP ones, but not much. Thus, it is preferable to allow the generated plans to split when added to the master problem.
The last runs (see table 3) were made to check how many plans it is interesting to generate. In the runs 28 – 31 the master is solved as an IP problem, and in the last two runs, the generated plans are fixed and the best chosen.
Settings / ResultsRun / Plans / DDEM shift / IDEM shift / Term crit / Opt Master LP / Opt Master IP / With plan nr
Monsteras / Morrum / Varo
28 / 50 / 2 / 5 / 2, 1% / 841200 / 844999 / 28 / 42 / 20
29 / 50 / 0 / 6 / 2, 1% / 836828 / 840214 / 16 / 16 / 30
30 / 50 / 2 / 4 / 2, 1% / 844472 / 847228 / 38 / 13 / 19
31 / 50 / 0 / 7 / 2, 1% / 845913 / 849597 / 38 / 24 / 22
32 / 50 / 2 / 5 / 2, 1% / 841200 / 888270 / 12
33 / 50 / 0 / 6 / 2, 1% / 841289 / 892296 / 25
Table 3: Runs with 50 generated plans.
As seen in table 3, it is possible to reduce the LP further, but the last generated plans are never used in the IP solution. Therefore it is not certain that a better solution would be found if more plans were generated. The last runs also clearly show the advantage of solving the master as an IP instead of fixing the production plans.
More interesting results should be obtained when the master is solved together with the ship-route generation part as well as this production plan generation problem. When more routes are available it is possible to estimate the demand shifts better and in that way get better production plans.
Appendix 1: The execution file used in the runs.
# Reseting and reading in the data and the model file.
reset;
reset options;
model model.txt;
data data.dat;
data prodn.dat;
data ships.dat;
option cplex_options
'mipdisplay = 3'
'mipinterval = 10000'
'mipgap 1e-2'
'mipsolutions 2'
#'time = 1200'
'probe = 1'
#'timing 1'
;
option show_stats 0;
# Below we define the number of decimals a solution is represented in.
# Both in displaying the solution and permenent in the program structures.
option solution_round 6;
# The main loop.
let TERMINATE:= 0;
param shift_domestic := 0;
param shift_international := 6;
for {t in PERIODS} {
for {p in PRODS} {
let TOTAL_DEMAND[t,p] :=
(if t <= last(PERIODS) - shift_domestic then
sum {d in DDEMS} DEMAND_DDEMS[d,p,t+shift_domestic]) +
(if t <= last(PERIODS) - shift_international then
sum {c in CUSTS} DEMAND_CUST[c,p,t+shift_international]);
}
}
repeat {
# Solving MASTER.
printf"** BEGIN SOLVING: MASTER **\n";
solve MASTER;
printf"** END SOLVING: MASTER **\n";
printf"** OPTIMAL VALUE: %1f **\n", totalcost;
# Extracting the dual variables from the MASTER problem.
let{m in MILLS, a in ASSORTS, t in PERIODS} FCA_PRICES[m,a,t]:= Storage_Mill_Ass[m,a,t].dual;
let{m in MILLS, p in PRODS, t in PERIODS} FCP_PRICES[m,p,t]:= Storage_Mill_Prod[m,p,t].dual;
let{m in MILLS} CON_PRICES[m]:= ProdPlanConvexity[m].dual;
#Solving subproblem.
printf"** SOLVING: SUBPROBLEM **\n";
solve SUBPROBLEM;
printf"** THE TOTAL REDUCED COST: %4f **\n", total_reduced_cost;
if total_reduced_cost < -0.0001
then {
printf"** SENDING THE COLUMNS TO MASTER! **\n";
for{m in MILLS}
{
let NUMPRODNPLANS[m]:= NUMPRODNPLANS[m] + 1;
for{r in RECIPES[m], t in PERIODS}
let PLANRECEIPEUSAGE[m, NUMPRODNPLANS[m], r, t]:= Y[m, r, t];
let C_PRODPLAN[m, NUMPRODNPLANS[m]] :=
sum{r in RECIPES[m], t in PERIODS} C_RECIPE[m, r] * Y[m, r, t] +
sum{r1 in RECIPES[m], r2 in RECIPES[m], t in PERIODS}
C_RECIPECHANGE[m, r1, r2] * X[m, r1, r2, t];
}
}
else
let TERMINATE:= 1;
printf"** iteration: %4f **\n", NUMPRODNPLANS["VARO"]-5;
} until (TERMINATE = 1 or NUMPRODNPLANS["VARO"] >= 50);
solve MASTER;
display prodplanusage;
printf"** OPTIMAL LP VALUE: %1f **\n", totalcost;
printf"** SOLVING IP MASTER PROBLEM **\n";
/*
option cplex_options
'mipdisplay = 3'
'mipinterval = 100'
'mipgap 1e-8'
'probe = 1'
;
redeclare var prodplanusage{j in MILLS, PRODPLANS[j]} binary;
solve MASTER;
display prodplanusage;*/
param raknare;
let raknare := 6;
repeat {
printf"** Production plan number: %1f **\n", (raknare-5);
fix prodplanusage["VARO", raknare] := 1;
fix prodplanusage["MONSTERAS", raknare] := 1;
fix prodplanusage["MORRUM", raknare+1] := 1;
solve MASTER;
unfix prodplanusage["VARO", raknare];
unfix prodplanusage["MONSTERAS", raknare];
unfix prodplanusage["MORRUM", raknare+1];
printf"** OPTIMAL IP VALUE: %1f **\n", totalcost;
let raknare := raknare + 1;
} until (raknare > NUMPRODNPLANS["VARO"]);