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