Decision Optimization Covid-19 Prototype Two
After a first prototype showing how Decision Optimization could help with the transfer of sick people over regions, a second prototype is shown here to illustrate how DO can help with the mask supply chain.
After the urgency of managing the capacity of reanimations rooms, moving cases from highly affected regions to less affected regions, the peak is almost reached, and the need is to focus on the how the population will be allowed to go out again. It seems that one of the main hypotheses and support for this end of confinement will be the availability of masks for everyone. Importations from China are starting again, and some few local factories are ramping up their production. It seems possible to set up a supply chain for masks for everyone.
In this example notebook, we imagine a supply chain based on importations, production in some few factories, important warehouses and distribution center in each French region (department).
Again, this is just a prototype to illustrate how Decision Optimization could be used, without any pretention that the model or its outcome are realistic. It should be used as a didactical illustrative example.
The following elements are taken into account:
- all productions, importations, consumptions and stocks are considered over
NB_PERIODS, - target areas are the French regions (departments), with their geographical location, and their consumptions for masks are considered to increase up to a peak and then stabilize,
- a few big warehouses in some geographies will stock the strategic reserves, with a given total capacity and initial stock,
- some national factories have some capacity to produce per period, which will increase over time,
- some importation occurs, only once per period, delivered to a unique warehouse per period to choose.
The example notebook is available here, and is made of 4 parts.
- loading and displaying geographical data on a map,
- simulating the evolution of consumption and production over time,
- creating and solving an optimization model
- displaying solutions.
Loading the data
A few csv files are used for factories, warehouses and departments. Again Folium is used to show this data on a map. The nearest warehouse for each factory and department, is also stored as data to be used later in the optimization model.
Simulating the production and consumption evolutions.
The data contains some current initial production and consumption for factories and departments.
The production is assumed to grow over time up to a maximum. The production for the 3 factories will look like this:
On the other hand the consumption is assumed to be more complex, with a growth up to a peak and then a small decrease and stabilization. The total consumption for all departments will then look like this:
Decision Optimization model
The optimization model is a very standard commodity flow problem. There are two types of sources: importation and production, and one type of sinks, the consumption at departments.
Decision variables include:
- production variables for each factory and period, all production goes fully and only to the nearest warehouse,
- importation variables for each warehouse and period, in addition to boolean variables indicating whether importation occurs for each warehouse and period,
- transportation variables between warehouses and departments,
- stock variables and under safety level variables.
transportation_periods = list(range(NB_PERIODS))fac_production_vars = mdl.integer_var_matrix(facs, all_periods, lb=0, name="factory production")import_vars = mdl.integer_var_matrix(wars, transportation_periods, lb=0, name="importation")
do_import_vars = mdl.binary_var_matrix(wars, transportation_periods, name="do importation")move_vars = mdl.integer_var_cube(wars, deps, transportation_periods, lb=0, name="transportation")war_stock_vars = mdl.integer_var_matrix(wars, all_periods, lb=0, name="warehouse stock")
war_under_safety_vars = mdl.integer_var_matrix(wars, all_periods, name="warehouse under satefy")
dep_stock_vars = mdl.integer_var_matrix(deps, all_periods, lb=0, name="departement stock")
dep_under_safety_vars = mdl.integer_var_matrix(deps, all_periods, name="departement under safety")mdl.print_information()
Constraints include:
- initial state constraints for stocks,
- safety stock definition constraints
- importation variables constraints,
- stock maximum capacity for departments and warehouses,
- conservation constraints over time for departments and warehouses
# Initial state
mdl.add_constraints(war_stock_vars[w, 0] == warehouses['initial'][w] for w in wars)
mdl.add_constraints(dep_stock_vars[d, 0] == departements['INITIAL'][d] for d in deps)# departement safety stock
mdl.add_constraints(mdl.max(0, int(departements['CAPACITY'][d]/SAFETY_LEVEL) - dep_stock_vars[d, p]) == dep_under_safety_vars[d,p] for d in deps for p in all_periods)# warehouse safety stock
mdl.add_constraints(mdl.max(0, int(warehouses['capacity'][w]/SAFETY_LEVEL) - war_stock_vars[w, p]) == war_under_safety_vars[w, p] for w in wars for p in all_periods)# import vars to do import vars links
mdl.add_constraints((import_vars[w, p] >= 1) == do_import_vars[w, p] for w in wars for p in transportation_periods)# import vars to one warehouse max per period
mdl.add_constraints(mdl.sum(do_import_vars[w, p] for w in wars) <= 1 for p in transportation_periods)# import vars max per period
mdl.add_constraints(mdl.sum(import_vars[w, p] for w in wars) <= MAX_IMPORTATION for p in transportation_periods)# departement stock capacity
mdl.add_constraints(dep_stock_vars[d, p] <= departements['CAPACITY'][d] for d in deps for p in all_periods)# factory production capacity
mdl.add_constraints(fac_production_vars[f, p] <= int(fac_production_capacity['capacity'][(fac_production_capacity['factory']==f) & (fac_production_capacity['period']==p)]) for f in facs for p in all_periods)# conservation constraints at warehouses
mdl.add_constraints(war_stock_vars[w, t+1] == war_stock_vars[w, t]
- mdl.sum(move_vars[w, d, t] for d in deps)
+ mdl.sum(fac_production_vars[f, t] for f in facs if factories['warehouse'][f]==w)
+ import_vars[w, t]
for w in wars for t in transportation_periods)# conservation constraints at departements
mdl.add_constraints(dep_stock_vars[d, t+1] == dep_stock_vars[d, t]
+ mdl.sum(move_vars[w, d, t] for w in wars)
- int(dep_consumption['consumption'][(dep_consumption['departement']==d) & (dep_consumption['period']==t)])
for d in deps for t in transportation_periods)mdl.print_information()
KPIs include:
- total importation,
- total production,
- total warehouse and department stock,
- total under safety level for warehouses and departments,
- transportation cost
The objective is to minimize under safety KPIs and (with a much lower weight) the transportation and importation costs.
total_importation = mdl.sum( import_vars[w, t] for w in wars for t in transportation_periods)
mdl.add_kpi(total_importation, 'total_importation')total_production = mdl.sum( fac_production_vars[f, t] for f in facs for t in all_periods)
mdl.add_kpi(total_production, 'total_production')total_warehouse_stock = mdl.sum( war_stock_vars[w, t] for w in wars for t in all_periods)
mdl.add_kpi(total_warehouse_stock, 'total_warehouse_stock')total_departement_stock = mdl.sum( dep_stock_vars[d, t] for d in deps for t in all_periods)
mdl.add_kpi(total_departement_stock, 'total_departement_stock')total_warehouse_under_safety = mdl.sum( war_under_safety_vars[w, t] for w in wars for t in all_periods)
mdl.add_kpi(total_warehouse_under_safety, 'total_warehouse_under_safety')total_departement_under_safety = mdl.sum( dep_under_safety_vars[d, t] for d in deps for t in all_periods)
mdl.add_kpi(total_departement_under_safety, 'total_departement_under_safety')transportation_cost = mdl.sum(d2w_distances[d][w] * move_vars[w, d, t] for w in wars for d in deps for t in transportation_periods)
mdl.add_kpi(transportation_cost, 'transportation_cost')mdl.minimize(10000*total_warehouse_under_safety + 10000*total_departement_under_safety + transportation_cost + total_importation)
mdl.print_information()
This simplistic model has around 20 000 variables and 20 000 constraints. This is still pretty small for a typical Decision Optimization model. With the CPLEX engine, problems are solve with millions of variables and constraints.
So the problem is solved pretty quickly.
CPXPARAM_Read_DataCheck 1
CPXPARAM_Threads 2
CPXPARAM_TimeLimit 20
Tried aggregator 2 times.
MIP Presolve eliminated 6650 rows and 2035 columns.
Aggregator did 6518 substitutions.
Reduced MIP has 3134 rows, 13452 columns, and 22291 nonzeros.
Reduced MIP has 1589 binaries, 10349 generals, 0 SOSs, and 3103 indicators.
Presolve time = 0.14 sec. (98.24 ticks)
Probing fixed 0 vars, tightened 3028 bounds.
Probing time = 0.04 sec. (10.01 ticks)
Tried aggregator 2 times.
MIP Presolve eliminated 0 rows and 2 columns.
Aggregator did 1 substitutions.
Reduced MIP has 3133 rows, 13449 columns, and 22288 nonzeros.
Reduced MIP has 1588 binaries, 10348 generals, 0 SOSs, and 3101 indicators.
Presolve time = 0.04 sec. (17.29 ticks)
Probing time = 0.01 sec. (3.71 ticks)
Clique table members: 15.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 2 threads.
Root relaxation solution time = 0.08 sec. (49.94 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
0 0 5.48064e+12 884 5.48064e+12 2424
0 0 5.48064e+12 474 Cuts: 883 2845
0 0 5.48064e+12 205 Cuts: 883 3151
0 0 5.48064e+12 60 Cuts: 620 3316
0 0 5.48064e+12 60 Flowpaths: 3 3320
* 0+ 0 5.48310e+12 5.48064e+12 0.04%
0 0 5.48064e+12 60 5.48310e+12 Flowpaths: 2 3322 0.04%
0 2 5.48064e+12 60 5.48310e+12 5.48064e+12 3322 0.04%
Elapsed time = 2.63 sec. (1282.59 ticks, tree = 0.02 MB, solutions = 1)
48 47 5.48064e+12 87 5.48310e+12 5.48064e+12 5455 0.04%
110 102 5.48073e+12 53 5.48310e+12 5.48064e+12 7503 0.04%
167 155 5.48239e+12 85 5.48310e+12 5.48064e+12 11793 0.04%
251 227 5.48141e+12 62 5.48310e+12 5.48064e+12 13649 0.04%
* 343+ 308 5.48224e+12 5.48064e+12 0.03%
347 332 5.48147e+12 28 5.48224e+12 5.48064e+12 14639 0.03%
391 370 5.48093e+12 80 5.48224e+12 5.48064e+12 17091 0.03%
* 395+ 369 5.48158e+12 5.48064e+12 0.02%
425 372 5.48066e+12 51 5.48158e+12 5.48064e+12 18788 0.02%
481 411 5.48139e+12 93 5.48158e+12 5.48064e+12 20961 0.02%
531 456 5.48087e+12 59 5.48158e+12 5.48064e+12 22020 0.02%
* 606+ 511 5.48156e+12 5.48064e+12 0.02%
* 646+ 556 5.48098e+12 5.48064e+12 0.01%
Cover cuts applied: 7
Implied bound cuts applied: 746
Flow path cuts applied: 11
Root node processing (before b&c):
Real time = 2.62 sec. (1280.66 ticks)
Parallel b&c, 2 threads:
Real time = 4.53 sec. (2800.44 ticks)
Sync time (average) = 0.32 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 7.16 sec. (4081.11 ticks)* KPI: total_importation = 150000000.000
* KPI: total_production = 20415342.000
* KPI: total_warehouse_stock = 1370635305.000
* KPI: total_departement_stock = 445073666.000
* KPI: total_warehouse_under_safety = 429364695.000
* KPI: total_departement_under_safety = 116926334.000
* KPI: transportation_cost = 17921196860.972
Display solution
The outcome of this model is the amount of importation per period and warehouse, the amount of production per factory and period, and the amount of transportation required between warehouses and departments.
To display the outcomes, Brunel visualizations is used to show the outcome as charts.
You can find the complete notebook here and the data here.
Conclusion
This new prototype illustrates yet another possible use of Decision Optimization to support government decisions in the Covid-19 crisis.
You can imagine many other problems which could benefit from Decision Optimization, such as the optimization of supply chain for tests.
