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