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

   file foliorep_sol.mos
   `````````````````````
   Modeling a MIP problem 
   to perform portfolio optimization.

   Same model as in foliomip3.mos.
   -- Infeasible model parameter values --
   -- Repairing infeasibilities --
   -- Storing solutions in an mpsol object --
   
  (c) 2014 Fair Isaac Corporation
      author: S.Heipcke, Sep. 2014, rev. Sep. 2022
*******************************************************!)

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.15                     ! Max. investment per ind. sector
  MAXVAL = 0.2                      ! Max. investment per share
  MINVAL = 0.1                      ! Min. investment per share
  MAXNUM = 4                        ! Max. number of different assets
  DATAFILE = "folio10.dat"          ! File with problem data
 end-parameters

 forward procedure print_sol
 forward procedure print_violated

 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
  buy: array(SHARES) of mpvar       ! 1 if asset is in portfolio, 0 otherwise
  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
  LimitNum: linctr                  ! Max. total number of assets
 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
 sum(s in SHARES) frac(s) = 1
 
! Upper bounds on the investment per share
 forall(s in SHARES) frac(s) <= MAXVAL

! Limit the total number of assets
 LimitNum:= sum(s in SHARES) buy(s) <= MAXNUM

 forall(s in SHARES) do
  buy(s) is_binary                  ! Turn variables into binaries
  frac(s) <= MAXVAL*buy(s)                 ! Linking the variables
  frac(s) >= MINVAL*buy(s)                 ! Linking the variables
 end-do


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

! Solve the problem
! maximize(XPRS_LIN,Return)
 maximize(Return)
 
!**** Infeasibility repair + reporting ****

 declarations
  Alrp,Agrp: array(linctr) of real  ! Selectors for LEG / GEQ constraints
  Albp,Aubp: array(mpvar) of real   ! Selector for lower / upper bounds on vars
  rstat: array(range) of string     ! Status message text
 end-declarations

 probstat:= getprobstat
 case probstat of
  XPRS_OPT: writeln("Problem solved")
  XPRS_INF: do 
             ! Must use the detailed infeasibility repair method since
	     ! only some constraints of each type may be relaxed
              Alrp(LimitRisk):=1
              forall(r in REGIONS) Alrp(LimitMaxReg(r)):=1
	      forall(t in TYPES) Alrp(LimitSec(t)):=1
              Alrp(LimitNum):=1
              forall(r in REGIONS) Agrp(LimitMinReg(r)):=1

	      rstat::(0..5)["relaxed optimum found", 
	        "relaxed problem infeasible", "relaxed problem unbounded", 
                "solution nonoptimal for original objective",
                "error", "numerical instability"]

              delta:=0.001
	      while (delta<10) do
	       setparam("XPRS_TIMELIMIT", 10)
!	       setparam("XPRS_VERBOSE", false)

               ! Option 'r': repairinfeas reports the relaxed constraints/bounds
                repairinfeas(Alrp, Agrp, Albp, Aubp, 'r', delta, "")
		probstat:= getprobstat
                write("delta = ", delta, ": ")
                case probstat of
                  XPRS_OPT: writeln("Relaxed problem solved")
                  XPRS_INF: writeln("Relaxed problem infeasible")
                  XPRS_OTH: writeln("Relaxed problem unbounded") 
                  XPRS_UNF: writeln("Repairinfeas has been terminated")
                  else writeln("Unknown problem status")
                end-case

               ! Display the relaxed solution
                print_sol
		
               ! Programmatic checking of constraint violations
                print_violated
		
                delta:= delta*10
	      end-do 
            end-do
  XPRS_OTH: writeln("Problem unbounded") 
  XPRS_UNF: writeln("Optimization unfinished")
  else writeln("Unknown problem status")
 end-case

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

!**** Display violated  constraints ****
 procedure print_violated
  declarations
   Solution: mpsol
   AllCtr: set of linctr
  end-declarations
 
  EPS:=1e-6
  
  writeln(" Violated (relaxed) constraints:")
  getloadedlinctrs(AllCtr)
  savesol(Solution)
  
  forall(c in AllCtr) 
   case gettype(c) of
    CT_EQ:  if abs(getsol(Solution,c))>EPS then
             writeln("   = constraint ", getname(c), " by ", getsol(Solution,c))
            end-if 
    CT_GEQ: if getsol(Solution,c)<-EPS then
             writeln("  >= constraint ", getname(c), " by ", getsol(Solution,c))
            end-if 
    CT_LEQ: if getsol(Solution,c)>EPS then
             writeln("  <= constraint ", getname(c), " by ", getsol(Solution,c))
            end-if 
   end-case
   writeln
 end-procedure

end-model