Dantzig-Wolfe decomposition: combining sequential and parallel solving

Dantzig-Wolfe decomposition (see [Teb01] for further detail) is a solution method for problems where, if a relatively small number of constraints were removed, the problem would fall apart into a number of independent problems. This means, it is possible to re-order the rows and columns of the constraint matrix as shown in Figure Coeffcient matrix with primal block angular structure, where non-zero coefficients only occur within the gray shaded areas. Such a primal block angular structure may become immediately apparent by visualizing a problem matrix. However, in most cases it will be necessary to re-organize the constraint definitions, grouping them by common index (sub)sets such as time periods, products, plant locations, and so on.

Figure 7: Coeffcient matrix with primal block angular structure

The constraints (including the objective function) involving variables of several or all subproblems are referred to as global constraints (also: common, linking, or master constraints). These constraints are used to form the master problem. The individual blocks on the diagonal of the coefficient matrix are solved as pricing subproblems, coordinated by the master problem. By solving the master problem we obtain a solution to the original problem. Since the master problem has a large number of variables (defined by the set of basic feasible solutions and unbounded directions of the pricing problems), we work with a restricted master problem over a small subset of the variables. The variables to enter the active set of variables of the restricted master problem are determined by solving the pricing subproblems. With objective functions for the pricing problems that are based on the dual values of the restricted master problem we can obtain that the objective function value at each extreme point is the reduced cost (or price) of the master problem variable associated with the extreme point.

For maximization problems, solving the modified pricing problems generates basic feasible solutions of maximum reduced cost. If the objective value at an extreme point is positive, then the associated master problem variable is added to the master problem; if the minimum objective value over all extreme points is negative, then no master problem variables exists to improve the current master problem solution.

The computational advantage of Dantzig-Wolfe decomposition arises from performing a significant amount of the computational work in the pricing problems that are roughly an order of magnitude smaller than the original problem and thus easier to solve. An aspect of the method that is of interest in the context of this paper is that the subproblems are independent of each other, and may therefore be solved concurrently.

A potential drawback of the decomposition approach is the huge size of the master problem, it has many more variables—though fewer constraints—than the original problem. In general it is not required to generate all variables explicitly but since the feasible region of the master problem is more complex than that of the original problem, the solution path may be longer. Furthermore, numerical problems may occur through the dynamic generation of variables of the master problem.

Many factors may influence the performance of a decomposition approach, so for a particular application computational experiments will be required to find out whether this solution method is suitable. Such tests may include different ways of decomposing a given problem. As a general rule for the definition of a problem decomposition, one should aim for few global constraints since it is important to be able to (re)solve the master problem quickly. In addition, the pricing problems should be constructed in such a way that they are well formed problems in their own right to avoid computational problems due to degeneracy.

Example problem: multi-item, multi-period production planning

The company Coco has two plants that can produce two types of cocoa powder. The first factory has a total capacity of 400 tons per month and the second of 500 tons per month. The marketing department has provided estimations for the maximum amount of every product that may be sold in each of the next four months, and also the expected sales prices. Several raw materials are used in the production with known raw material requirements per ton of finished products. Finished products and raw material may be stored at the factories from one tim period to the next, incurring a given cost per ton and per time period. The raw material storage capacity at the factories is limited to 300 tons. How should the two plants be operated during the planning period to maximize the total profit?

Original model

Let PRODS be the set of finished products, FACT the set of factories, RAW the set of raw materials, and TIMES={1,...,NT} the set of time periods under consideration.

We define decision variables make_pft representing the quantity of product p made at factory f during time period t. Furthermore, to model the transition from one time period to the next and to account for the different types of cost incurred, we need several other sets of variables: sell_pft, the amount of product p sold at factory f in period t; buy_rft the amount of raw material r bought at f in t; and finally pstock_pft and rstock_rft (both defined for t=1,...,NT+1) respectively the amount of product and raw material held in stock at factory f at the beginning of period t.

Let further MXSELL_pt be the maximum sales quantity of product p in period t, MXMAKE_f the capacity limit of factory f, and MXRSTOCK the raw material storage capacity.

Let IPSTOCK_pf be the quantity of product p initially held in stock at factory f and IRSTOCK_rf the initial stock level of raw material r. We denote by CPSTOCK and CRSTOCK respectively the unit cost for storing finished product and raw material.

The objective function of maximizing the total profit is to maximize the sales revenues (REV_p), minus the cost of production (CMAKE_pf), buying raw material (CBUY_rt), and storing finished products and raw material.

The possibility to store products between time periods gives rise to three sets of constraints: the inventory balance constraints for finished products (PBal_pft) and for raw material (RBal_rft), and the limit on the raw material storage capacity (MxRStock_ft). The stock pstock_p,f,t+1 of product p at the beginning of t+1 is given by the stock at the beginning of t plus the production in t reduced by the amount sold on t. The stock rstock_r,f,t+1 of raw material r at the beginning of t+1 is given by the corresponding stock at the beginning of t plus the amount bought in t reduced by the quantity used in production during t.

∀ p ∈ PRODS, ∀ f ∈ FACT, ∀ t ∈ TIME: PBal_pft:= pstock_p,f,t+1 = pstock_pft + make_pft - sell_pft
∀ r ∈ RAW, ∀ f ∈ FACT, ∀ t ∈ TIME:
RBal_rft:= rstock_r,f,t+1 = rstock_rft + buy_rft -

∑

_{p ∈ PRODS}

REQ_pr · make_pft
∀ f ∈ FACT, ∀ t ∈ {2,...,NT+1}: MxRStock_ft:=

∑

_{r ∈ RAW}

rstock_rft ≤ MXRSTOCK

We further have two sets of capacity constraints: the production capacity of the factories is limited (constraints MxMake_ft) and we can only sell up to a given maximum amount of finished products per time period (constraints MxSell_ft).

Below the complete mathematical model is given. The initial stock levels (t=1) of finished products and raw material are fixed to the given values.

maximize

∑

_{p ∈ PRODS}

∑

_{f ∈ FACT}

∑

_{t ∈ TIME}

REV_pt · sell_pft
-

∑

_{p ∈ PRODS}

∑

_{f ∈ FACT}

∑

_{t ∈ TIME}

CMAKE_pf · make_pft -

∑

_{r ∈ RAW}

∑

_{f ∈ FACT}

∑

_{t ∈ TIME}

CBUY_rt · buy_rft
-

∑

_{p ∈ PRODS}

∑

_{f ∈ FACT}

^NT+1

∑

_{t =2}

CPSTOCK · pstock_pft -

∑

_{r ∈ RAW}

∑

_{f ∈ FACT}

^NT+1

∑

_{t =2}

CRSTOCK · rstock_rft
∀ p ∈ PRODS, ∀ f ∈ FACT, ∀ t ∈ TIME: PBal_pft:= pstock_p,f,t+1 = pstock_pft + make_pft - sell_pft
∀ r ∈ RAW, ∀ f ∈ FACT, ∀ t ∈ TIME:
RBal_rft:= rstock_r,f,t+1 = rstock_rft + buy_rft -

∑

_{p ∈ PRODS}

REQ_pr · make_pft
∀ p ∈ PRODS, ∀ t ∈ TIME: MxSell_pt:=

∑

_{f ∈ FACT}

sell_pft ≤ MXSELL_p
∀ f ∈ FACT, ∀ t ∈ TIME: MxMake_ft:=

∑

_{p ∈ PRODS}

make_pft ≤ MXMAKE_f
∀ f ∈ FACT, ∀ t ∈ {2,...,NT+1}: MxRStock_ft:=

∑

_{r ∈ RAW}

rstock_rft ≤ MXRSTOCK
∀ p ∈ PRODS, ∀ f ∈ FACT: pstock_pf1 = IPSTOCK_pf
∀ r ∈ RAW, ∀ f ∈ FACT: rstock_rf1 = IRSTOCK_rf
∀ p ∈ PRODS, ∀ f ∈ FACT, ∀ t ∈ TIME: make_pft≥ 0, sell_pft≥ 0
∀ r ∈ RAW, ∀ f ∈ FACT, ∀ t ∈ TIME: buy_rft ≥ 0
∀ p ∈ PRODS, ∀ f ∈ FACT, ∀ t ∈ {1,...,NT+1}: pstock_pft ≥ 0
∀ r ∈ RAW, ∀ f ∈ FACT, ∀ t ∈ {1,...,NT+1}: rstock_rft ≥ 0

Problem decomposition

We now decompose the problem described above according to production locations. Notice that this is not the only way of decomposing this problem: we may just as well choose a decomposition by products or by time periods. However, both of these choices result in a larger number of global constraints than the decomposition by factories, meaning that the master problem may become more difficult to solve.

For every factory f we obtain the following subproblem (including the sales limit constraints MxSell in the form of bounds in the submodels is not required, but may lead to better, that is, more exact or faster solutions).

maximize

∑

_{p ∈ PRODS}

∑

_{t ∈ TIME}

REV_pt · sell_pt
-

∑

_{p ∈ PRODS}

∑

_{t ∈ TIME}

CMAKE_p · make_pt -

∑

_{r ∈ RAW}

∑

_{t ∈ TIME}

CBUY_rt · buy_rt
-

∑

_{p ∈ PRODS}

^NT+1

∑

_{t =2}

CPSTOCK · pstock_pt -

∑

_{r ∈ RAW}

^NT+1

∑

_{t =2}

CRSTOCK · rstock_rt
∀ p ∈ PRODS, ∀ t ∈ TIME: PBal_pt:= pstock_p,t+1 = pstock_pt + make_pt - sell_pt
∀ r ∈ RAW, ∀ t ∈ TIME: RBal_rt:= rstock_r,t+1 = rstock_rt + buy_rt -

∑

_{p ∈ PRODS}

REQ_pr · make_pt
∀ t ∈ TIME: MxMake_t:=

∑

_{p ∈ PRODS}

make_pt ≤ MXMAKE
∀ t ∈ {2,...,NT+1}: MxRStock_t:=

∑

_{r ∈ RAW}

rstock_rt ≤ MXRSTOCK
∀ p ∈ PRODS, ∀ t ∈ TIME: MxSell_pt:= sell_pt ≤ MXSELL_p
∀ p ∈ PRODS: pstock_p1 = IPSTOCK_p
∀ r ∈ RAW: rstock_r1 = IRSTOCK_r
∀ p ∈ PRODS, ∀ t ∈ TIME: make_pt≥ 0, sell_pt≥ 0
∀ r ∈ RAW, ∀ t ∈ TIME: buy_rt ≥ 0
∀ p ∈ PRODS, ∀ t ∈ {1,...,NT+1}: pstock_pt ≥ 0
∀ r ∈ RAW, ∀ t ∈ {1,...,NT+1}: rstock_rt ≥ 0

The master problem only contains a single set of global constraints, namely the sales limit constraints MxSell (for clarity's sake, the non-negativity constraints are left out here).

In the decomposition algorithm, the decision variables in the master problem are expressed via the solutions (proposals) generated by the subproblems, such as

∀ p ∈ PRODS, ∀ f ∈ FACT, ∀ t ∈ TIME: sell_pft =

^nPROP_f

∑

_k=1

Prop_sell_pftk·weight_fk

where Prop_sell_pftk stands for the solution value of variable sell_pt in the k^th proposal generated by subproblem f and Prop_cost_fk is the objective value this proposal. For every subproblem f we need to add a convexity constraint Convex_f on the weight_fk variables.

maximize

∑

_{f ∈ FACT}

^nPROP_f

∑

_k=1

Prop_cost_fk·weight_fk
∀ p ∈ PRODS, ∀ t ∈ TIME: MxSell_pt:=

∑

_{f ∈ FACT}

^nPROP_f

∑

_k=1

Prop_sell_pftk·weight_fk ≤ MXSELL_p
∀ f ∈ FACT: Convex_f:=

^nPROP_f

∑

_k=1

weight_fk = 1
∀ f ∈ FACT, ∀ k ∈ 1,...,nPROP_f: weight_fk≥ 0

We shall refer to this problem as the modified master problem. Without going any further into technical detail we simply remark that a correspondence exists between the solution of the original problem and those of the modified master problem.

Implementation

The decomposition algorithm has several phases:

Phase 1: generation of a set of proposals to obtain a feasible solution to the modified master problem.
Phase 2: optimization of the modified master problem.
Phase 3: calculate the solution to the original problem.

The suproblems solved in phases 1 and 2 take as their objective functions sums of the dual values from the modfied master problem. To avoid starting off with an empty master problem and hence random dual information that may be misleading we add a crash Phase 0 to our implementation that generates one proposal for each subproblem with the original objective function.

Master model

Below follows the body of the master model file. The definitions of the function bodies will be shown later in Section Master problem subroutines.

model "Coco3 Master"
 uses "mmxprs","mmjobs","mmsystem"

 parameters
  DATAFILE = "coco2.dat"
  ALG = 0                               ! 0: stop phase with 1st failed subpb.
                                        ! 1: stop when all subprob.s fail
 end-parameters

 forward procedure process_sub_result
 forward procedure solve_master(phase:integer)
 forward procedure process_master_result
 forward function calc_solution:real
 forward procedure print_solution

 declarations
  PHASE_0=2                           ! Event codes sent to submodels
  PHASE_1=3
  PHASE_2=4
  PHASE_3=5
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_FAILED=7
  EVENT_READY=8
  NPROD, NFACT, NRAW, NT: integer
 end-declarations

 initializations from DATAFILE
  NPROD NFACT NRAW NT
 end-initializations

 declarations
  PRODS = 1..NPROD                    ! Range of products (p)
  FACT = 1..NFACT                     !          factories (f)
  RAW = 1..NRAW                       !          raw materials (r)
  TIME = 1..NT                        !          time periods (t)

  nIter: integer                      ! Iteration counter
  nPROP: array(FACT) of integer       ! Counters of proposals from subprob.s
 end-declarations

!**** Master problem ****
 declarations
  MXSELL: array(PRODS,TIME) of real   ! Max. amount of p that can be sold
  excessS: mpvar                      ! Violation of sales/buying limits
  weight: array(FACT,range) of mpvar  ! weights for proposals
  MxSell: array(PRODS,TIME) of linctr ! Sales limit constraints
  Convex: array(FACT) of linctr       ! Convexity constraints
  Price_convex: array(FACT) of real     ! Dual price on convexity constraints
  Price_sell: array(PRODS,TIME) of real ! Dual price on sales limits
 end-declarations

 initializations from DATAFILE
  MXSELL
 end-initializations

!**** Submodels ****
 declarations
  submod: array(FACT) of Model        ! One subproblem per factory
  Stopped: set of integer
 end-declarations

 res:= compile("g","cocoSubF.mos")    ! Compile the submodel file
 forall(f in FACT) do                 ! Load & run one submodel per product
  Price_convex(f):= 1
  load(submod(f), "cocoSubF.bim")
  submod(f).uid:= f
  run(submod(f), "Factory=" + f + ",DATAFILE=" + DATAFILE)
  wait                                ! Wait for child model to be ready
  ev:=getnextevent
  if ev.class=EVENT_END then
   writeln("Cannot start all necessary models - aborting")
   exit(1)
  end-if
 end-do

!**** Phase 0: Crash ****
 nIter:=1; finished:=false
 writeln("\nPHASE 0 -- Iteration ", nIter); fflush

 forall(f in FACT)                    ! Start solving all submodels (Phase 1)
  send(submod(f), PHASE_0, 0)

 forall(f in FACT) do
  wait                                ! Wait for child (termination) events
  ev:= getnextevent
  if getclass(ev)=EVENT_SOLVED then
   process_sub_result                 ! Add new proposal to master problem
  elif getclass(ev)=EVENT_FAILED then
   finished:= true
  end-if
 end-do

 if finished then
  writeln("Problem is infeasible")
  exit(1)
 end-if

 solve_master(1)                      ! Solve the updated Ph. 1 master problem
 process_master_result                ! Store initial pricing data for submodels


!**** Phase 1: proposal generation (feasibility) ****
 repeat
  noimprove:= 0
  nIter+=1
  writeln("\nPHASE 1 -- Iteration ", nIter); fflush

  forall(f in FACT)                   ! Start solving all submodels (Phase 1)
   send(submod(f), PHASE_1, Price_convex(f))

  forall(f in FACT) do
   wait                               ! Wait for child (termination) events
   ev:= getnextevent
   if getclass(ev)=EVENT_SOLVED then
    process_sub_result                ! Add new proposal to master problem
   elif getclass(ev)=EVENT_FAILED then
    noimprove += 1
   end-if
  end-do

  if noimprove = NFACT then
   writeln("Problem is infeasible")
   exit(2)
  end-if
  if ALG=0 and noimprove > 0 then
   writeln("No improvement by some subproblem(s)")
   break
  end-if

  solve_master(1)                     ! Solve the updated Ph. 1 master problem
  if getobjval>0.00001 then
   process_master_result              ! Store new pricing data for submodels
  end-if
 until getobjval <= 0.00001


!**** Phase 2: proposal generation (optimization) ****
 writeln("\n**** PHASE 2 ****")
 finished:=false
 repeat
  solve_master(2)                     ! Solve Phase 2 master problem
  process_master_result               ! Store new pricing data for submodels

  nIter+=1
  writeln("\nPHASE 2 -- Iteration ", nIter); fflush

  forall(f in FACT)                   ! Start solving all submodels (Phase 2)
   send(submod(f), PHASE_2, Price_convex(f))

  forall(f in FACT) do
   wait                               ! Wait for child (termination) events
   ev:= getnextevent
   if getclass(ev)=EVENT_SOLVED then
    process_sub_result                ! Add new proposal to master problem
   elif getclass(ev)=EVENT_FAILED then
    if ALG=0 then
     finished:=true                   ! 1st submodel w/o prop. stops phase 2
    else
     Stopped += {ev.fromuid}          ! Stop phase 2 only when no submodel
                                      ! generates a new proposal
    end-if
   end-if
  end-do

  if getsize(Stopped) = NFACT then finished:= true; end-if
 until finished

 solve_master(2)                      ! Re-solve master to integrate
                                      ! proposal(s) from last ph. 2 iteration

!**** Phase 3: solution to the original problem ****
 writeln("\n**** PHASE 3 ****")
 forall(f in FACT) do
  send(submod(f), PHASE_3, 0)         ! Stop all submodels
  wait
  dropnextevent
 end-do

 writeln("Total Profit=", calc_solution)
 print_solution
end-model

The initial declarations block of this model defines a certain number of event codes that are used to identify the messages sent between this master model and its child (sub)models. The same declarations need to be repeated in the child models.

Single factory model

The model file cocoSubF.mos with the definition of the subproblems has the following contents.

model "Coco Subproblem (factory based decomp.)"
 uses "mmxprs", "mmjobs"

 parameters
  Factory = 0
  TOL = 0.00001
  DATAFILE = "coco3.dat"
 end-parameters

 forward procedure process_solution

 declarations
  PHASE_0=2                             ! Event codes sent to submodels
  PHASE_1=3
  PHASE_2=4
  PHASE_3=5
  EVENT_SOLVED=6                        ! Event codes sent by submodels
  EVENT_FAILED=7
  EVENT_READY=8
  NPROD, NFACT, NRAW, NT: integer
 end-declarations

 send(EVENT_READY,0)                    ! Model is ready (= running)

 initializations from DATAFILE
  NPROD NFACT NRAW NT
 end-initializations

 declarations
  PRODS = 1..NPROD                      ! Range of products (p)
  FACT = 1..NFACT                       !          factories (f)
  RAW = 1..NRAW                         !          raw materials (r)
  TIME = 1..NT                          !          time periods (t)

  REV: array(PRODS,TIME) of real        ! Unit selling price of products
  CMAKE: array(PRODS,FACT) of real      ! Unit cost to make product p
                                        ! at factory f
  CBUY: array(RAW,TIME) of real         ! Unit cost to buy raw materials
  REQ: array(PRODS,RAW) of real         ! Requirement by unit of product p
                                        ! for raw material r
  MXSELL: array(PRODS,TIME) of real     ! Max. amount of p that can be sold
  MXMAKE: array(FACT) of real           ! Max. amount factory f can make
                                        ! over all products
  IPSTOCK: array(PRODS,FACT) of real    ! Initial product stock levels
  IRSTOCK: array(RAW,FACT) of real      ! Initial raw material stock levels
  CPSTOCK: real                         ! Unit cost to store any product p
  CRSTOCK: real                         ! Unit cost to store any raw mat. r
  MXRSTOCK: real                        ! Raw material storage capacity

  make: array(PRODS,TIME) of mpvar      ! Amount of products made at factory
  sell: array(PRODS,TIME) of mpvar      ! Amount of product sold from factory
  buy: array(RAW,TIME) of mpvar         ! Amount of raw material bought
  pstock: array(PRODS,1..NT+1) of mpvar ! Product stock levels at start
                                        ! of period t
  rstock: array(RAW,1..NT+1) of mpvar   ! Raw material stock levels
                                        ! at start of period t

  sol_make: array(PRODS,TIME) of real   ! Amount of products made
  sol_sell: array(PRODS,TIME) of real   ! Amount of product sold
  sol_buy: array(RAW,TIME) of real      ! Amount of raw mat. bought
  sol_pstock: array(PRODS,1..NT+1) of real ! Product stock levels
  sol_rstock: array(RAW,1..NT+1) of real   ! Raw mat. stock levels

  Profit: linctr                        ! Profit of proposal
  Price_sell: array(PRODS,TIME) of real ! Dual price on sales limits
 end-declarations

 initializations from DATAFILE
  CMAKE REV CBUY REQ MXSELL MXMAKE
  IPSTOCK IRSTOCK MXRSTOCK CPSTOCK CRSTOCK
 end-initializations

! Product stock balance
 forall(p in PRODS,t in TIME)
  PBal(p,t):= pstock(p,t+1) = pstock(p,t) + make(p,t) - sell(p,t)

! Raw material stock balance
 forall(r in RAW,t in TIME)
  RBal(r,t):= rstock(r,t+1) =
   rstock(r,t) + buy(r,t) - sum(p in PRODS) REQ(p,r)*make(p,t)

! Capacity limit
 forall(t in TIME)
  MxMake(t):= sum(p in PRODS) make(p,t) <= MXMAKE(Factory)

! Limit on the amount of prod. p to be sold
 forall(p in PRODS,t in TIME) sell(p,t) <= MXSELL(p,t)

! Raw material stock limit
 forall(t in 2..NT+1)
  MxRStock(t):= sum(r in RAW) rstock(r,t) <= MXRSTOCK

! Initial product and raw material stock levels
 forall(p in PRODS) pstock(p,1) = IPSTOCK(p,Factory)
 forall(r in RAW) rstock(r,1) = IRSTOCK(r,Factory)

! Total profit
 Profit:=
  sum(p in PRODS,t in TIME) REV(p,t) * sell(p,t) -           ! revenue
  sum(p in PRODS,t in TIME) CMAKE(p,Factory) * make(p,t) -   ! prod. cost
  sum(r in RAW,t in TIME) CBUY(r,t) * buy(r,t) -             ! raw mat.
  sum(p in PRODS,t in 2..NT+1) CPSTOCK * pstock(p,t) -       ! p storage
  sum(r in RAW,t in 2..NT+1) CRSTOCK * rstock(r,t)           ! r storage

! (Re)solve this model until it is stopped by event "PHASE_3"
 repeat
  wait
  ev:= getnextevent
  Phase:= getclass(ev)
  if Phase=PHASE_3 then               ! Stop the execution of this model
   break
  end-if
  Price_convex:= getvalue(ev)         ! Get new pricing data

  if Phase<>PHASE_0 then
   initializations from "bin:shmem:pricedata"
    Price_sell
   end-initializations
  end-if

 ! (Re)solve this model
  if Phase=PHASE_0 then
   maximize(Profit)
  elif Phase=PHASE_1 then
   maximize(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) + Price_convex)
  else        ! PHASE 2
   maximize(
    Profit - sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) -
     Price_convex)
  end-if

  writeln("Factory ", Factory, " - Obj: ", getobjval,
          " Profit: ", getsol(Profit), " Price_sell: ",
          getsol(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) ),
          " Price_convex: ", Price_convex)
  fflush

  if getobjval > TOL then             ! Solution found: send values to master
   process_solution
  elif getobjval <= TOL then          ! Problem is infeasible (Phase 0/1) or
   send(EVENT_FAILED,0)               ! no improved solution found (Phase 2)
  else
   send(EVENT_READY,0)
  end-if
 until false

!-----------------------------------------------------------
! Process solution data
 procedure process_solution
  forall(p in PRODS,t in TIME) do
   sol_make(p,t):= getsol(make(p,t))
   sol_sell(p,t):= getsol(sell(p,t))
  end-do
  forall(r in RAW,t in TIME) sol_buy(r,t):= getsol(buy(r,t))
  forall(p in PRODS,t in 1..NT+1) sol_pstock(p,t):= getsol(pstock(p,t))
  forall(r in RAW,t in 1..NT+1) sol_rstock(r,t):= getsol(rstock(r,t))
  Prop_cost:= getsol(Profit)
  send(EVENT_SOLVED,0)

  initializations to "mempipe:noindex,sol"
   Factory
   sol_make sol_sell sol_buy sol_pstock sol_rstock
   Prop_cost
  end-initializations
 end-procedure
end-model

The child models are re-solved until they receive the PHASE_3 code. At every iteration they write their solution values to memory so that these can be processed by the master model.

Master problem subroutines

The following three subroutines of the master model recover the solutions produced by the subproblems (process_sub_result), re-solve the master problem (solve_master), and communicate the solution of the master problem to its child models (process_master_result).

 declarations
  Prop_make: array(PRODS,FACT,TIME,range) of real ! Amount of products made
  Prop_sell: array(PRODS,FACT,TIME,range) of real ! Amount of product sold
  Prop_buy: array(RAW,FACT,TIME,range) of real    ! Amount of raw mat. bought
  Prop_pstock: array(PRODS,FACT,1..NT+1,range) of real ! Product stock levels
  Prop_rstock: array(RAW,FACT,1..NT+1,range) of real   ! Raw mat. stock levels
  Prop_cost: array(FACT,range) of real  ! Cost/profit of each proposal
 end-declarations

 procedure process_sub_result
  declarations
   f: integer                         ! Factory index
                                      ! Solution values of the proposal:
   sol_make: array(PRODS,TIME) of real      ! Amount of products made
   sol_sell: array(PRODS,TIME) of real      ! Amount of product sold
   sol_buy: array(RAW,TIME) of real         ! Amount of raw mat. bought
   sol_pstock: array(PRODS,1..NT+1) of real ! Product stock levels
   sol_rstock: array(RAW,1..NT+1) of real   ! Raw mat. stock levels
   pc: real                           ! Cost of the proposal
  end-declarations

 ! Read proposal data from memory
  initializations from "mempipe:noindex,sol"
   f
   sol_make sol_sell sol_buy sol_pstock sol_rstock
   pc
  end-initializations

 ! Add the new proposal to the master problem
  nPROP(f)+=1
  create(weight(f,nPROP(f)))
  forall(p in PRODS,t in TIME) do
   Prop_make(p,f,t,nPROP(f)):= sol_make(p,t)
   Prop_sell(p,f,t,nPROP(f)):= sol_sell(p,t)
  end-do
  forall(r in RAW,t in TIME) Prop_buy(r,f,t,nPROP(f)):= sol_buy(r,t)
  forall(p in PRODS,t in 1..NT+1) Prop_pstock(p,f,t,nPROP(f)):= sol_pstock(p,t)
  forall(r in RAW,t in 1..NT+1) Prop_rstock(r,f,t,nPROP(f)):= sol_rstock(r,t)
  Prop_cost(f,nPROP(f)):= pc
  writeln("Sol. for factory ", f, ":\n  make:   ", sol_make, "\n  sell:   ",
           sol_sell, "\n  buy:    ", sol_buy, "\n  pstock: ", sol_pstock,
	   "\n  rstock: ", sol_rstock)
 end-procedure

!-----------------------------------------------------------
 procedure solve_master(phase: integer)
  forall(f in FACT)
   Convex(f):= sum (k in 1..nPROP(f)) weight(f,k) = 1

  if phase=1 then
   forall(p in PRODS,t in TIME)
    MxSell(p,t):=
     sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) -
      excessS <= MXSELL(p,t)
   minimize(excessS)
  else
   forall(p in PRODS,t in TIME)
    MxSell(p,t):=
     sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) <=
      MXSELL(p,t)
   maximize(sum(f in FACT, k in 1..nPROP(f)) Prop_cost(f,k) * weight(f,k))
  end-if

  writeln("Master problem objective: ", getobjval)
  write("  Weights:")
  forall(f in FACT,k in 1..nPROP(f)) write(" ", getsol(weight(f,k)))
  writeln
 end-procedure

!-----------------------------------------------------------
 procedure process_master_result
  forall(p in PRODS,t in TIME) Price_sell(p,t):=getdual(MxSell(p,t))
  forall(f in FACT) Price_convex(f):=getdual(Convex(f))

  initializations to "bin:shmem:pricedata"
   Price_sell
  end-initializations
 end-procedure

Finally, the master model is completed by two subroutines for calculating the solution to the original problem (calc_solution), Phase 3 of the decomposition algorithm, and printing out the solution (print_solution). The solution to the original problem is obtained from the solution values of the modfied master problem and the proposals generated by the subproblems.

 declarations
  REV: array(PRODS,TIME) of real      ! Unit selling price of products
  CMAKE: array(PRODS,FACT) of real    ! Unit cost to make product p
                                      ! at factory f
  CBUY: array(RAW,TIME) of real       ! Unit cost to buy raw materials
  COPEN: array(FACT) of real          ! Fixed cost of factory f being
                                      ! open for one period
  CPSTOCK: real                       ! Unit cost to store any product p
  CRSTOCK: real                       ! Unit cost to store any raw mat. r

  Sol_make: array(PRODS,FACT,TIME) of real ! Solution value (products made)
  Sol_sell: array(PRODS,FACT,TIME) of real ! Solution value (product sold)
  Sol_buy: array(RAW,FACT,TIME) of real    ! Solution value (raw mat. bought)
  Sol_pstock: array(PRODS,FACT,1..NT+1) of real ! Sol. value (prod. stock)
  Sol_rstock: array(RAW,FACT,1..NT+1) of real   ! Sol. value (raw mat. stock)
 end-declarations

 initializations from DATAFILE
  CMAKE REV CBUY CPSTOCK CRSTOCK COPEN
 end-initializations

 function calc_solution: real
  forall(p in PRODS,f in FACT,t in TIME) do
   Sol_sell(p,f,t):=
    sum(k in 1..nPROP(f)) Prop_sell(p,f,t,k) * getsol(weight(f,k))
   Sol_make(p,f,t):=
    sum(k in 1..nPROP(f)) Prop_make(p,f,t,k) * getsol(weight(f,k))
  end-do
  forall(r in RAW,f in FACT,t in TIME) Sol_buy(r,f,t):=
    sum(k in 1..nPROP(f)) Prop_buy(r,f,t,k) * getsol(weight(f,k))
  forall(p in PRODS,f in FACT,t in 1..NT+1) Sol_pstock(p,f,t):=
   sum(k in 1..nPROP(f)) Prop_pstock(p,f,t,k) * getsol(weight(f,k))
  forall(r in RAW,f in FACT,t in 1..NT+1) Sol_rstock(r,f,t):=
   sum(k in 1..nPROP(f)) Prop_rstock(r,f,t,k) * getsol(weight(f,k))

  returned:=
   sum(p in PRODS,f in FACT,t in TIME) REV(p,t) * Sol_sell(p,f,t) -
   sum(p in PRODS,f in FACT,t in TIME) CMAKE(p,f) * Sol_make(p,f,t) -
   sum(r in RAW,f in FACT,t in TIME) CBUY(r,t) * Sol_buy(r,f,t) -
   sum(p in PRODS,f in FACT,t in 2..NT+1) CPSTOCK * Sol_pstock(p,f,t) -
   sum(r in RAW,f in FACT,t in 2..NT+1) CRSTOCK * Sol_rstock(r,f,t)
 end-function

!-----------------------------------------------------------
 procedure print_solution
  writeln("Finished products:")
  forall(f in FACT) do
   writeln("Factory ", f, ":")
   forall(p in PRODS) do
    write("  ", p, ":    ")
    forall(t in TIME) write(strfmt(Sol_make(p,f,t),6,1), "(",
                            strfmt(Sol_pstock(p,f,t+1),5,1), ")")
    writeln
   end-do
  end-do

  writeln("Raw material:")
  forall(f in FACT) do
   writeln("Factory ", f, ":")
   forall(r in RAW) do
    write("  ", r, ": ")
    forall(t in TIME) write(strfmt(Sol_buy(r,f,t),6,1), "(",
                            strfmt(Sol_rstock(r,f,t+1),5,1), ")")
    writeln
   end-do
  end-do

  writeln("Sales:")
  forall(f in FACT) do
   writeln("Factory ", f, ":")
   forall(p in PRODS) do
    write("  ", p, ": ")
    forall(t in TIME) write(strfmt(Sol_sell(p,f,t),4))
    writeln
   end-do
  end-do

  writeln("\nComputation time: ", gettime)
 end-procedure

Results

For the test cases that have been tried the solutions produced by our decomposition algorithm are close to the optimal solution, but the latter is not always reached. The reason behind this is that the decomposition algorithm is a sequence of iterations that may accumulate errors at different levels---lowering the tolerances used as stopping criterion in the submodels most of the time does not improve the solution. However, the configuration of the decomposition algorithm itself shows some impact on the solution: in phases 1 and 2 one may choose, for instance, to stop once the first submodel returns no improvement or continue until no more proposals are generated. Generating more proposals sometimes helps finding a better solution, but it also increases the number of times (sub)problems are solved and hence prolongates the solving time.

Contents

Index

Glossary

Search Results