Alternative implementation with multiple problems

Since our decomposition algorithm is formed by a series of sequential optimization runs it is possible to implement the whole approach as a single model defining several problems (instead of several model files, each with a single problem). This alternative implementation requires significantly less effort for data handling (problems simply share data) and coordination of solving (problem solving in a single model always is sequential). As a result, the total Mosel model code is shorter. However, by collecting all problems into a single model, the testing of (sub)problems in isolation is made somewhat more complicated—in the previous version we can simply run one of the (sub)models as a stand-alone model, a feature that certainly is helpful when developing large applications.

Master model

In the master model, all code related to (sub)model handling is replaced by the somewhat simpler handling of (sub)problems that are defined and solved within the same model as the master problem. Here we just display the part of the declarations that are new or modified, followed by the setup of the subproblems and the (unchanged) solution algorithm with the definition of the `start solve' subroutines.

For every submodel we now have two new subroutines, define_*prob and solve_*prob that replace the separate model files. These subroutines are lsited in the following sections.

 declarations
  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

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

! **** 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 

Submodel 1: fixed continuous variables

The bounds and integrality conditions for the integer subproblem are stated once (define_intprob). Each time the integer subproblem is solved (solve_intprob) a new constraint gets added. The solving routine also stores the solution values and the problem status into global data structures.

! 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

Submodel 2: fixed integer variables

The setup routine for the continuous problem (define_contprob) only defines the objective function. The constraints are re-defined each time the problem gets solved (solve_contprob).

! 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

Submodel 0: start solution

The setup routine for the dual problem (define_dualprob) states the constraints for this model. The solving subroutines solve_dualprob implements two cases, depending on whether we use this model with a dummy objective for finding a start solution (Step 0) or in the second part of the solution algorithm where the objective is defined with the current solution values (Step 2).

! 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

!-----------------------------------------------------------
! (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 

Similarly to the previous implementation as a separate model, we have moved the storing of the solution values into a subroutine (save_dualsolution). Notice that it is not necessary to switch problems using with prob do in this subroutine because it gets called from within a subproblem.

! Process dual solution data
 procedure save_dualsolution(prob:mpproblem)
 ! 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