(!******************************************************
   Mosel Example Problems
   ======================

   file folioinfeas.mos
   ````````````````````
   Modeling an LP problem 
   to perform portfolio optimization.

   Linear version of the model in foliomip3.mos.
   -- Infeasible model parameter values --
   -- Handling infeasibility through auxiliary variables --
   
  (c) 2008 Fair Isaac Corporation
      author: S.Heipcke, Dec. 2008, rev. Sep. 2018
*******************************************************!)

model "Portfolio optimization with MIP"
 uses "mmxprs"

 parameters
  MAXRISK = 1/3                     ! Max. investment into high-risk values
  MINREG = 0.3                      ! Min. investment per geogr. region
  MAXREG = 0.5                      ! Max. investment per geogr. region
  MAXSEC = 0.1                      ! Max. investment per ind. sector
  MAXVAL = 0.2                      ! Max. investment per share
  MINVAL = 0.1                      ! Min. investment per share
  DATAFILE = "folio250.dat"         ! File with problem data
 end-parameters

 forward procedure print_sol

 declarations
  SHARES: set of string              ! Set of shares
  RISK: set of string                ! Set of high-risk values among shares
  REGIONS: set of string             ! Geographical regions
  TYPES: set of string               ! Share types (ind. sectors)
  LOC: array(REGIONS) of set of string ! Sets of shares per geogr. region
  RET: array(SHARES) of real         ! Estimated return in investment
  SEC: array(TYPES) of set of string ! Sets of shares per industry sector
 end-declarations

 initializations from DATAFILE
  RISK RET LOC SEC
 end-initializations

 public declarations
  frac: array(SHARES) of mpvar      ! Fraction of capital used per share
  Return: linctr                    ! Total return
  LimitRisk: linctr                 ! Max. percentage of high-risk values
  LimitMinReg,LimitMaxReg: array(REGIONS) of linctr  ! Min/max perc. per region
  LimitSec: array(TYPES) of linctr  ! Max. percentage per industry sector
  TotalOne: linctr                  ! Spend all the capital
 end-declarations

! Objective: total return
 Return:= sum(s in SHARES) RET(s)*frac(s) 

! Limit the percentage of high-risk values
 LimitRisk:= sum(s in RISK) frac(s) <= MAXRISK

! Limits on geographical distribution
 forall(r in REGIONS) do
  LimitMinReg(r):= sum(s in LOC(r)) frac(s) >= MINREG
  LimitMaxReg(r):= sum(s in LOC(r)) frac(s) <= MAXREG
 end-do 

! Diversification across industry sectors
 forall(t in TYPES) LimitSec(t):= sum(s in SEC(t)) frac(s) <= MAXSEC

! Spend all the capital
 TotalOne:= sum(s in SHARES) frac(s) = 1
 
! Upper bounds on the investment per share
 forall(s in SHARES) frac(s) <= MAXVAL

! Display Optimizer log
 setparam("XPRS_verbose", true)

! Solve the problem (using primal Simplex to minimize infeasibilities)
 maximize(XPRS_PRI,Return)

 declarations
  mx,mdj: set of mpvar
  mslack,mdual: set of linctr
  devRisk: mpvar
  devMinReg,devMaxReg: array(REGIONS) of mpvar
  devSec: array(TYPES) of mpvar
 end-declarations 

 if getprobstat=XPRS_INF then
 ! Print some details
  getinfeas(mx, mslack, mdual, mdj)
  writeln("Infeasible variables:")
  forall(v in mx) writeln("  ",getname(v))
  forall(v in mdj) writeln("  ",getname(v))
  writeln("Infeasible constraints:")
  forall(c in mslack) writeln("  ",getname(c))
  forall(c in mdual) writeln("  ",getname(c))
 
  writeln("#### Original problem infeasible. Adding deviation variables ####")

 ! Add deviation variables to make problem solvable
  LimitRisk -= devRisk
  forall(r in REGIONS) do
   LimitMinReg(r) += devMinReg(r)
   LimitMaxReg(r) -= devMaxReg(r)
   devMaxReg(r) <= MAXREG/2  ! Only allow small deviations
  end-do
  forall(t in TYPES) do
   LimitSec(t) -= devSec(t)
   devSec(t)<= MAXSEC/2              ! Only allow small deviations
  end-do 
 
  Penalty:= devRisk + sum(r in REGIONS) (devMinReg(r)+devMaxReg(r)) +
            sum(t in TYPES) devSec(t)
 
 ! Resolve the problem with penalty terms added to the objective
  maximize(Return - 10*Penalty)

  if getprobstat=XPRS_INF then
   writeln("No solution after relaxation")
   exit(1)
  end-if
  
  writeln("Constraint violations:")
  writeln("  Constraint            Activity  Deviation  Bound(s)")
  writeln("  Risk                  ", 
          strfmt(getsol(sum(s in RISK) frac(s)),-10), 
          strfmt(devRisk.sol,-10), " (", MAXRISK, ")")
  forall(r in REGIONS)
   writeln("  Region ", strfmt(r,-11), "    ", 
           strfmt(getsol(sum(s in LOC(r)) frac(s)),-10), 
           strfmt(devMaxReg(r).sol-devMinReg(r).sol,-10), 
           " (", MINREG, "-",MAXREG,")")
  forall(t in TYPES)
   writeln("  Sector ", strfmt(t,-14), " ", 
           strfmt(getsol(sum(s in SEC(t)) frac(s)),-10),  
           strfmt(devSec(t).sol,-10), " (", MAXSEC, ")")    
 end-if

 print_sol

! Solution printing
 procedure print_sol
  writeln("Solution ", getparam("XPRS_MIPSOLS"))
  writeln("Total return: ", getsol(Return))
  forall(s in SHARES | getsol(frac(s))>0)
   writeln(s, ": ", getsol(frac(s))*100, "%")
 end-procedure
end-model