Initializing help system before first use

Benders decomposition: sequential solving of several different submodels


Type: Programming
Rating: 5 (difficult)
Description: Benders decomposition is a method for solving large MIP problems. The model implementation shows the following features:
  • iterative sequence of concurrent solving of a set of subproblems,
  • data exchange between several models via shared memory, and
  • coordination of several models via events.
An implementation using a single model is also presented (benders_single.mos).
File(s): benders_main.mos, benders_dual.mos, benders_cont.mos, benders_int.mos, benders_single.mos
Data file(s): bprob12.dat, bprob33.dat


benders_main.mos
(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_main.mos
   `````````````````````
   Benders decomposition for solving a simple MIP.
   - Master model - 

   *** ATTENTION: This model will return an error if ***
   *** no more than one Xpress licence is available. ***

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2006, rev. Dec. 2017
*******************************************************!)
        
model "Benders (master model)"
 uses "mmxprs", "mmjobs", "mmsystem"

 parameters
  NCTVAR = 3  ! 1
  NINTVAR = 3 ! 2
  NC = 4      ! 3
  BIGM = 1000
  ALG = 1                             ! 1: Use Benders decomposition (dual)
                                      ! 2: Use Benders decomposition (primal)
  DATAFILE = "bprob33.dat"  ! "bprob12.dat"
 end-parameters

 forward procedure start_solution
 forward procedure solve_primal_int(ct: integer)
 forward procedure solve_cont
 forward function eval_solution: boolean
 forward procedure print_solution

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_1=3
  STEP_2=4
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_INFEAS=7
  EVENT_READY=8

  CtVars = 1..NCTVAR                  ! Continuous variables
  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  A: array(Ctrs,CtVars) of integer    ! Coeff.s of continuous variables
  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  C: array(CtVars) of integer         ! Obj. coeff.s of continuous variables
  D: array(IntVars) of integer        ! Obj. coeff.s of discrete variables
  Ctr: array(Ctrs) of linctr          ! Constraints of orig. problem
  CtrD: array(CtVars) of linctr       ! Constraints of dual problem
  MC: array(range) of linctr          ! Constraints generated by alg.

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_x: array(CtVars) of real        ! Solution of primal prob. (cont.)
  sol_y: array(IntVars) of real       ! Solution of primal prob. (integers)
  sol_obj: real                       ! Objective function value (primal)

  RM: range                           ! Model indices 
  stepmod: dynamic array(RM) of Model ! Submodels
 end-declarations

 initializations from DATAFILE
  A B b C D
 end-initializations

! **** Submodels ****

 initializations to "bin:shmem:probdata"   ! Save data for submodels
  A  B  b  C  D  
 end-initializations

! Compile + load all submodels
 if compile("benders_int.mos")<>0 then exit(1); end-if
 create(stepmod(1)); load(stepmod(1), "benders_int.bim")
 if compile("benders_dual.mos")<>0 then exit(2); end-if

 if ALG=1 then
  create(stepmod(2)); load(stepmod(2), "benders_dual.bim")
 else
  create(stepmod(0)); load(stepmod(0), "benders_dual.bim")
  if compile("benders_cont.mos")<>0 then exit(3); end-if
  create(stepmod(2)); load(stepmod(2), "benders_cont.bim")
  run(stepmod(0), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
 end-if
                                      ! Start the execution of the submodels
 run(stepmod(1), "NINTVAR=" + NINTVAR + ",NC=" + NC)
 run(stepmod(2), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)

 forall(m in RM) do
  wait                                ! Wait for "Ready" messages
  ev:= getnextevent
  if getclass(ev) <> EVENT_READY then
   writeln("Error occurred in a subproblem")
   exit(4)
  end-if
 end-do

! **** Solution algorithm ****

 start_solution                       ! 0. Initial solution for cont. var.s
 ct:= 1
 repeat
 writeln("\n**** Iteration: ", ct)
  solve_primal_int(ct)                ! 1. Solve problem with fixed cont.  
  solve_cont                          ! 2. Solve problem with fixed int.
  ct+=1
 until eval_solution                  !    Test for optimality 
 print_solution                       ! 3. Retrieve and display the solution


! **** Cleaning up temporary files ****
 fdelete("benders_int.bim")
 fdelete("benders_dual.bim")
 if ALG<>1 then fdelete("benders_cont.bim"); end-if
 fdelete("shmem:probdata")
 fdelete("shmem:sol")
 
!----------------------------------------------------------- 
! Produce an initial solution for the dual problem using a dummy objective
 procedure start_solution
  if ALG=1 then                       ! Start the problem solving
   send(stepmod(2), STEP_0, 0)
  else
   send(stepmod(0), STEP_0, 0)
  end-if
  wait                                ! Wait for the solution
  ev:=getnextevent
  if getclass(ev)=EVENT_INFEAS then
   writeln("Problem is infeasible")
   exit(6) 
  end-if
 end-procedure

!----------------------------------------------------------- 
! Solve a modified version of the primal problem, replacing continuous 
! variables by the solution of the dual
 procedure solve_primal_int(ct: integer)  
  send(stepmod(1), STEP_1, ct)        ! Start the problem solving
  wait                                ! Wait for the solution
  ev:=getnextevent
  sol_obj:= getvalue(ev)              ! Store objective function value

  initializations from "bin:shmem:sol"  ! Retrieve the solution
   sol_y 
  end-initializations
 end-procedure

!----------------------------------------------------------- 
! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of y
 procedure solve_cont
  send(stepmod(2), STEP_2, 0)         ! Start the problem solving
  wait                                ! Wait for the solution
  dropnextevent

  initializations from "bin:shmem:sol"  ! Retrieve the solution
   sol_u 
  end-initializations
 end-procedure

!----------------------------------------------------------- 
 function eval_solution: boolean
  write("Test optimality: ", sol_obj - sum(i in IntVars) D(i)*sol_y(i),
          " = ", sum(j in Ctrs) sol_u(j)* (b(j) - 
	                            sum(i in IntVars) B(j,i)*sol_y(i)) )
  returned:= ( sol_obj - sum(i in IntVars) D(i)*sol_y(i) = 
      sum(j in Ctrs) sol_u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i)) )
  writeln(if(returned, " : true", " : false"))
 end-function

!----------------------------------------------------------- 
 procedure print_solution  
  ! Retrieve results
  initializations from "bin:shmem:sol"
   sol_x 
  end-initializations

  forall(m in RM) stop(stepmod(m))     ! Stop all submodels
   
  write("\n**** Solution (Benders): ", sol_obj, "\n  x: ")
  forall(i in CtVars) write(sol_x(i), " ")
  write("  y: ")
  forall(i in IntVars) write(sol_y(i), " ")
  writeln
 end-procedure

end-model

benders_dual.mos
(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_dual.mos
   `````````````````````
   Benders decomposition for solving a simple MIP.
   - Dual problem in the continuous variables
     for start solution and Step 2 of the algorithm -

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2006, rev. Jan. 2013
*******************************************************!)
        
model "Benders (dual problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3
  NINTVAR = 3
  NC = 4 
  BIGM = 1000
 end-parameters

 forward procedure save_solution  

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_2=4
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_INFEAS=7
  EVENT_READY=8

  CtVars = 1..NCTVAR                  ! Continuous variables
  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  A: array(Ctrs,CtVars) of integer    ! Coeff.s of continuous variables
  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  C: array(CtVars) of integer         ! Obj. coeff.s of continuous variables

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_x: array(CtVars) of real        ! Solution of primal prob. (cont.)
  sol_y: array(IntVars) of real       ! Solution of primal prob. (integers)
 
  u: array(Ctrs) of mpvar             ! Dual variables
 end-declarations

 initializations from "bin:shmem:probdata"
  A  B  b  C
 end-initializations
  
 forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)

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

! (Re)solve this model until it is stopped by event "STEP_3"
 repeat
  wait
  ev:= getnextevent
  Alg:= getclass(ev)

  if Alg=STEP_0 then                  ! Produce an initial solution for the 
                                      ! dual problem using a dummy objective
   maximize(sum(j in Ctrs) u(j))
   if(getprobstat = XPRS_INF) then 
    writeln("Problem is infeasible")
    send(EVENT_INFEAS,0)              ! Problem is infeasible
   else
    write("**** Start solution: ")
    save_solution
   end-if

  else                                ! STEP 2: Solve the dual problem for 
                                      ! given solution values of y
   initializations from "bin:shmem:sol"
    sol_y 
   end-initializations

   Obj:= sum(j in Ctrs) u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i))
   maximize(XPRS_PRI, Obj)

   if(getprobstat=XPRS_UNB) then
    BigM:= sum(j in Ctrs) u(j) <= BIGM
    maximize(XPRS_PRI, Obj)
   end-if
   
   write("Step 2: ")
   save_solution                      ! Write solution to memory
   BigM:= 0                           ! Reset the 'BigM' constraint
  end-if

 until false

!-----------------------------------------------------------
! Process solution data
 procedure save_solution  
 ! Store values of u and x
  forall(j in Ctrs) sol_u(j):= getsol(u(j))
  forall(i in CtVars) sol_x(i):= getdual(CtrD(i))

  initializations to "bin:shmem:sol"
   sol_u  sol_x 
  end-initializations
  
  send(EVENT_SOLVED, getobjval)

  write(getobjval, "\n  u: ")
  forall(j in Ctrs) write(sol_u(j), " ")
  write("\n  x: ")
  forall(i in CtVars) write(getdual(CtrD(i)), " ")
  writeln
  fflush
 end-procedure
   
end-model

benders_cont.mos
(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_cont.mos
   `````````````````````
   Benders decomposition for solving a simple MIP.
   - Solve the primal problem for given solution values 
     of the integer variable (Step 2 of the algorithm) -

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2006, rev. Jan. 2013
*******************************************************!)
        
model "Benders (continuous problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3
  NINTVAR = 3
  NC = 4 
  BIGM = 1000
 end-parameters

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_2=4
  STEP_3=5
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_READY=8

  CtVars = 1..NCTVAR                  ! Continuous variables
  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  A: array(Ctrs,CtVars) of integer    ! Coeff.s of continuous variables
  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  C: array(CtVars) of integer         ! Obj. coeff.s of continuous variables
  Ctr: array(Ctrs) of linctr          ! Constraints of orig. problem

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_x: array(CtVars) of real        ! Solution of primal prob. (cont.)
  sol_y: array(IntVars) of real       ! Solution of primal prob. (integers)
 
  x: array(CtVars) of mpvar           ! Continuous variables
 end-declarations

 initializations from "bin:shmem:probdata"
  A  B  b  C  
 end-initializations
  
 Obj:= sum(i in CtVars) C(i)*x(i)

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

! (Re)solve this model until it is stopped by event "STEP_3"
 repeat
  wait
  dropnextevent

  initializations from "bin:shmem:sol"
   sol_y 
  end-initializations

  forall(j in Ctrs)
   Ctr(j):= sum(i in CtVars) A(j,i)*x(i) + 
            sum(i in IntVars) B(j,i)*sol_y(i) >= b(j) 

  minimize(Obj)                       ! Solve the problem

 ! Store values of u and x
  forall(j in Ctrs) sol_u(j):= getdual(Ctr(j))
  forall(i in CtVars) sol_x(i):= getsol(x(i))

  initializations to "bin:shmem:sol"
   sol_u  sol_x
  end-initializations

  send(EVENT_SOLVED, getobjval)

  write("Step 2: ", getobjval, "\n  u: ")
  forall(j in Ctrs) write(sol_u(j), " ")
  write("\n  x: ")
  forall(i in CtVars) write(getsol(x(i)), " ")
  writeln
  fflush

 until false
   
end-model

benders_int.mos
(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_int.mos
   ````````````````````
   Benders decomposition for solving a simple MIP.
   - Solve the primal problem for given solution values 
     of the continuous variables (Step 1 of the algorithm) -

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2006, rev. Jan. 2013
*******************************************************!)
        
model "Benders (integer problem)"
 uses "mmxprs", "mmjobs"

 parameters
  NINTVAR = 3
  NC = 4 
  BIGM = 1000
 end-parameters

 declarations
  STEP_0=2                            ! Event codes sent to submodels
  STEP_1=3
  EVENT_SOLVED=6                      ! Event codes sent by submodels
  EVENT_READY=8

  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  D: array(IntVars) of integer        ! Obj. coeff.s of discrete variables
  MC: array(range) of linctr          ! Constraints generated by alg.

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_y: array(IntVars) of real       ! Solution of primal prob.
 
  y: array(IntVars) of mpvar          ! Discrete variables
  z: mpvar                            ! Objective function variable 
 end-declarations

 initializations from "bin:shmem:probdata"
  B  b  D
 end-initializations

 z is_free                            ! Objective is a free variable  
 forall(i in IntVars) y(i) is_integer ! Integrality condition
 forall(i in IntVars) y(i) <= BIGM    ! Set (large) upper bound

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

 repeat
  wait
  ev:= getnextevent
  ct:= integer(getvalue(ev))

  initializations from "bin:shmem:sol"
   sol_u
  end-initializations
    
  ! Add a new constraint
  MC(ct):= z >= sum(i in IntVars) D(i)*y(i) +  
                sum(j in Ctrs) sol_u(j)*(b(j) - sum(i in IntVars) B(j,i)*y(i))
    
  minimize(z)
  
  ! Store solution values of y
  forall(i in IntVars) sol_y(i):= getsol(y(i))

  initializations to "bin:shmem:sol"
   sol_y
  end-initializations
  
  send(EVENT_SOLVED, getobjval)

  write("Step 1: ", getobjval, "\n  y: ")
  forall(i in IntVars) write(sol_y(i), " ")
  
  write("\n  Slack: ")
  forall(j in 1..ct) write(getslack(MC(j)), " ")
  writeln
  fflush

 until false
   
end-model

benders_single.mos
(!*******************************************************
   Mosel Example Problems 
   ======================

   file benders_single.mos
   ```````````````````````
   Benders decomposition for solving a simple MIP.
   - Single model containing several subproblems - 

   (c) 2009 Fair Isaac Corporation
       author: S. Heipcke, Jan. 2009, rev. Dec. 2017
*******************************************************!)
        
model "Benders (single model)"
 uses "mmxprs", "mmjobs"

 parameters
  NCTVAR = 3  ! 1
  NINTVAR = 3 ! 2
  NC = 4      ! 3
  BIGM = 1000
  ALG = 1                             ! 1: Use Benders decomposition (dual)
                                      ! 2: Use Benders decomposition (primal)
  DATAFILE = "bprob33.dat"  ! "bprob12.dat"
 end-parameters

 forward procedure start_solution
 forward procedure solve_primal_int(ct: integer)
 forward procedure solve_cont
 forward function eval_solution: boolean
 forward procedure print_solution

 forward procedure define_intprob(prob:mpproblem)
 forward function solve_intprob(prob:mpproblem, ct:integer): real
 forward procedure define_dualprob(prob:mpproblem)
 forward function solve_dualprob(prob:mpproblem, Alg:integer): real
 forward procedure define_contprob(prob:mpproblem)
 forward function solve_contprob(prob:mpproblem): real
 
 declarations
  STEP_0=2                            ! Codes sent to subproblems
  STEP_1=3
  STEP_2=4
  STAT_SOLVED=6                       ! Status codes returned by subproblems
  STAT_INFEAS=7
  STAT_READY=8

  CtVars = 1..NCTVAR                  ! Continuous variables
  IntVars = 1..NINTVAR                ! Discrete variables
  Ctrs = 1..NC                        ! Set of constraints (orig. problem)

  A: array(Ctrs,CtVars) of integer    ! Coeff.s of continuous variables
  B: array(Ctrs,IntVars) of integer   ! Coeff.s of discrete variables
  b: array(Ctrs) of integer           ! RHS values
  C: array(CtVars) of integer         ! Obj. coeff.s of continuous variables
  D: array(IntVars) of integer        ! Obj. coeff.s of discrete variables
  Ctr: array(Ctrs) of linctr          ! Constraints of orig. problem
  CtrD: array(CtVars) of linctr       ! Constraints of dual problem
  MC: array(range) of linctr          ! Constraints generated by alg.
  Obj: linctr                         ! Continuous objective

  y: array(IntVars) of mpvar          ! Discrete variables
  z: mpvar                            ! Objective function variable 
  u: array(Ctrs) of mpvar             ! Dual variables
  x: array(CtVars) of mpvar           ! Continuous variables

  sol_u: array(Ctrs) of real          ! Solution of dual problem
  sol_x: array(CtVars) of real        ! Solution of primal prob. (cont.)
  sol_y: array(IntVars) of real       ! Solution of primal prob. (integers)
  sol_obj: real                       ! Objective function value (primal)

  RM: range                           ! Problem indices 
  stepprob: dynamic array(RM) of mpproblem    ! Subproblems
  status: dynamic array(mpproblem) of integer ! Subproblem status
 end-declarations

 initializations from DATAFILE
  A B b C D
 end-initializations

! **** Subproblems ****
                                      ! Create and define submodels
 create(stepprob(1)); define_intprob(stepprob(1))
 
 if ALG=1 then
  create(stepprob(2)); define_dualprob(stepprob(2))
 else
  create(stepprob(0)); define_dualprob(stepprob(0))
  create(stepprob(2)); define_contprob(stepprob(2))
 end-if

! **** Solution algorithm ****

 start_solution                       ! 0. Initial solution for cont. var.s
 ct:= 1
 repeat
 writeln("\n**** Iteration: ", ct)
  solve_primal_int(ct)                ! 1. Solve problem with fixed cont.  
  solve_cont                          ! 2. Solve problem with fixed int.
  ct+=1
 until eval_solution                  !    Test for optimality 
 print_solution                       ! 3. Retrieve and display the solution

!----------------------------------------------------------- 
! Produce an initial solution for the dual problem using a dummy objective
 procedure start_solution
  num:= if(ALG=1, 2, 0)
  res:=solve_dualprob(stepprob(num), STEP_0)  ! Start the problem solving

  if status(stepprob(num))=STAT_INFEAS then
   writeln("Problem is infeasible")
   exit(6) 
  end-if
 end-procedure

!----------------------------------------------------------- 
! Solve a modified version of the primal problem, replacing continuous 
! variables by the solution of the dual
 procedure solve_primal_int(ct: integer)  
  sol_obj:= solve_intprob(stepprob(1), ct)
 end-procedure

!----------------------------------------------------------- 
! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of y
 procedure solve_cont
  if ALG=1 then                       ! Start the problem solving
   res:= solve_dualprob(stepprob(2), STEP_2)
  else 
   res:= solve_contprob(stepprob(2))
  end-if
 end-procedure

!----------------------------------------------------------- 
 function eval_solution: boolean
  write("Test optimality: ", sol_obj - sum(i in IntVars) D(i)*sol_y(i),
          " = ", sum(j in Ctrs) sol_u(j)* (b(j) - 
	                            sum(i in IntVars) B(j,i)*sol_y(i)) )
  returned:= ( sol_obj - sum(i in IntVars) D(i)*sol_y(i) = 
      sum(j in Ctrs) sol_u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i)) )
  writeln(if(returned, " : true", " : false"))
 end-function

!----------------------------------------------------------- 
 procedure print_solution  
      ! TODO: Reset all submodels
   
  write("\n**** Solution (Benders): ", sol_obj, "\n  x: ")
  forall(i in CtVars) write(sol_x(i), " ")
  write("  y: ")
  forall(i in IntVars) write(sol_y(i), " ")
  writeln
 end-procedure

!----------------------------------------------------------- 
! Define the integer problem 
 procedure define_intprob(prob:mpproblem)
  with prob do
   z is_free                            ! Objective is a free variable  
   forall(i in IntVars) y(i) is_integer ! Integrality condition
   forall(i in IntVars) y(i) <= BIGM    ! Set (large) upper bound
  end-do
  status(prob):= STAT_READY
 end-procedure 
 
! Solve the integer problem 
 function solve_intprob(prob:mpproblem, ct:integer): real
  with prob do
   status(prob):= STAT_READY
    
   ! Add a new constraint
   MC(ct):= z >= sum(i in IntVars) D(i)*y(i) +  
                 sum(j in Ctrs) sol_u(j)*(b(j) - sum(i in IntVars) B(j,i)*y(i))
    
   minimize(z)
  
   ! Store solution values of y
   forall(i in IntVars) sol_y(i):= getsol(y(i))
  
   returned:=getobjval
   status(prob):= STAT_SOLVED

   write("Step 1: ", getobjval, "\n  y: ")
   forall(i in IntVars) write(sol_y(i), " ")  
   write("\n  Slack: ")
   forall(j in 1..ct) write(getslack(MC(j)), " ")
   writeln
   fflush

  end-do
 end-function
 
!----------------------------------------------------------- 
! Define the dual problem 
 procedure define_dualprob(prob:mpproblem)
  with prob do
   forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)
  end-do
  status(prob):= STAT_READY
 end-procedure 

! Process dual solution data
 procedure save_dualsolution(prob:mpproblem)
 ! There is no need for 'with mpproblem do' here since this subroutine
 ! is called from within a subproblem
 
 ! Store values of u and x
  forall(j in Ctrs) sol_u(j):= getsol(u(j))
  forall(i in CtVars) sol_x(i):= getdual(CtrD(i))
  
  status(prob):= STAT_SOLVED

  write(getobjval, "\n  u: ")
  forall(j in Ctrs) write(sol_u(j), " ")
  write("\n  x: ")
  forall(i in CtVars) write(getdual(CtrD(i)), " ")
  writeln
  fflush
 end-procedure

! (Re)solve the dual problem
 function solve_dualprob(prob:mpproblem, Alg:integer): real
  with prob do
   status(prob):= STAT_READY

   if Alg=STEP_0 then                 ! Produce an initial solution for the 
                                      ! dual problem using a dummy objective
    maximize(sum(j in Ctrs) u(j))
    if(getprobstat = XPRS_INF) then 
     writeln("Problem is infeasible")
     status(prob):= STAT_INFEAS       ! Problem is infeasible
    else
     write("**** Start solution: ")
     save_dualsolution(prob)
     returned:= getobjval
    end-if

   else                               ! STEP 2: Solve the dual problem for 
                                      ! given solution values of y
    Obj:= sum(j in Ctrs) u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i))
    maximize(XPRS_PRI, Obj)

    if(getprobstat=XPRS_UNB) then
     BigM:= sum(j in Ctrs) u(j) <= BIGM
     maximize(XPRS_PRI, Obj)
    end-if
   
    write("Step 2: ")
    save_dualsolution(prob)           ! Save solution values
    returned:= getobjval
    BigM:= 0                          ! Reset the 'BigM' constraint
   end-if

  end-do
 end-function

!-----------------------------------------------------------
! Define the continuous primal problem 
 procedure define_contprob(prob:mpproblem)
  with prob do
   Obj:= sum(i in CtVars) C(i)*x(i)
  end-do
  status(prob):= STAT_READY
 end-procedure 
 
! Solve the continuous problem 
 function solve_contprob(prob:mpproblem): real
  with prob do
   status(prob):= STAT_READY
 
  ! (Re)define and solve this problem
   forall(j in Ctrs)
    Ctr(j):= sum(i in CtVars) A(j,i)*x(i) + 
             sum(i in IntVars) B(j,i)*sol_y(i) >= b(j) 

   minimize(Obj)                       ! Solve the problem

  ! Store values of u and x
   forall(j in Ctrs) sol_u(j):= getdual(Ctr(j))
   forall(i in CtVars) sol_x(i):= getsol(x(i))

   returned:=getobjval
   status(prob):= STAT_SOLVED

   write("Step 2: ", getobjval, "\n  u: ")
   forall(j in Ctrs) write(sol_u(j), " ")
   write("\n  x: ")
   forall(i in CtVars) write(getsol(x(i)), " ")
   writeln
   fflush
  
  end-do
 end-function

!-----------------------------------------------------------

end-model