model DC_Location
 uses "mmxprs", "mmsheet", "mmnl"
 !setparam ("XPRS_VERBOSE", true)
 setparam("XPRS_MAXTIME",300)
 setparam("XPRS_MIPRELSTOP",0.01)
 	
 	
 declarations
 	
   ! --- Indices ---

D = 1..60
R = 1..60
I = 1..2

  ! --- Paramètre ---
  Fd    : real   
  T     : integer   
  Hr    : real  
  hd    : real  
  Kd    : integer    
  LoadingRate : real
  TotalLoad, TotalCapacity : real

  ! --- Distances & capacités ---
  Ddr   : array(R,D) of real    
  Cdr   : real  

  ! --- Demandes ---
  lambda : array(R,I) of real 	! lambda(r,i) unifié
  lambda_irSA1 : array(R) of real     
  lambda_irSA2 : array(R) of real
  
  Lambda_r : array(R) of real
  Lambda_i : array(I) of real
  
  ! --- Incertitude ---
  sigma  : array(R,I) of real    
  z_alpha : real   

  ! --- Lead times ---
  ! LTdr et LTdi sont égaux et lus directement depuis l'Excel
  LT  : real 

  ! --- Lot sizes ---
  Qdr   : array(D,R) of real     
  Qstar : array(D,R) of real      ! permet de calculer Cap%
  Lambda_ref : real
  CdrFixed : array(1..5,1..4) of real
  QlistCap : array(1..7200) of real
  CapPercent : integer
  IndexCap : integer
  k : integer
  temp : real
  nbContraints, nbNonContraints : integer


  ! --- Variables ---
  y    : array(D) of mpvar  
  x    : array(D,R) of mpvar  

  v1   : array(D,I) of mpvar  
  v2   : array(D,I) of mpvar
  
  
  ! --- Objectif ---
  OpenCost, FixedOpenCost, VariableOpenCost, FrequencyMeasure : linctr
  TransCost, StockRCost, StockDCCost, SafetyRCost, SafetyDCCost : linctr
  TotalCostOut, OpenCostOut, TransCostOut, StockRCostOut : real !export
  StockDCCostOut, SafetyRCostOut, SafetyDCCostOut : real
  
  ! --- Calibration de Cope ---
  TargetFixedShare : real
  FixedOpenCostOut, VariableOpenCostOut, FrequencyMeasureOut : real

  BaselineFixedSA1, BaselineFixedSA2 : real
  BaselineFrequencySA1, BaselineFrequencySA2 : real

  BaselineFixedTotal, BaselineFrequencyTotal : real
  COpeCalculated : real

  CalibrationMode : boolean
  
  ShareIndex : integer
  share : integer
  
  ObservedFixedShareOut : real
  
  
   ! --- Expériences ---
  MapID : integer
  NbDC : integer
  SCENARIO: integer
  
  ! --- Export Excel ---
  RowExport : integer
  CapIndex : integer
  FdIndex : integer
  
  !--- Lecture Maps ---
  MapFile : string
  
  
  !--- Contraintes ---
  C_v1     : array(D,I) of nlctr
  C_v2     : array(D,I) of nlctr
  C_assign : array(R) of linctr
  C_open   : array(D,R) of linctr
  
  !--- affichage GAP & TPS ---
  TimeOut, GapOut : real
  
  
! --- Frquence des expéditions ---
COpe : real
  
 end-declarations

procedure RunModel
	

   writeln("Run : Map=",MapID,
            " Cap=",CapPercent,
            " Fd=",Fd,
            " Scenario=",SCENARIO)
   MapFile := "Maps" + strfmt(MapID,0) + "B_60R_50_250_500.xlsx"
   writeln("Fichier lu : ", MapFile)
   
 initializations from "mmsheet.excel:noindex;" + MapFile
 	
 	!Fd as "[DonneesGenerales$J14]"
 	T as "[DonneesGenerales$N3]"
 	Hr as "[DonneesGenerales$J2]"
 	hd as "[DonneesGenerales$J6]"
 	Kd as "[DonneesGenerales$J10]"
 	
 	Ddr as "[Distances$B2:$BI61]"
 	!Cdr as "[DonneesGenerales$G6]"
 	
 	lambda_irSA1 as "[DonneesGenerales$C3:$C62]"
 	lambda_irSA2 as "[DonneesGenerales$D3:$D62]"
 	!lambda_irCOOP as "[DonneesGenerales$E3:$E62]"
 	
 	sigma as "[DonneesGenerales$M7:$N66]"
 	z_alpha as "[DonneesGenerales$G10]"
 	
 	! LTdr et LTdi sont égaux et lus directement depuis l'Excel
 	LT as "[DonneesGenerales$G3]"
 	
 	
 end-initializations

! ============================================================
! CHOIX DU SCENARIO
! ============================================================

forall(r in R) do
  if SCENARIO = 1 then
    lambda(r,1) := lambda_irSA1(r)
    lambda(r,2) := 0

  elif SCENARIO = 2 then
    lambda(r,1) := 0
    lambda(r,2) := lambda_irSA2(r)

  elif SCENARIO = 3 then
    lambda(r,1) := lambda_irSA1(r)
    lambda(r,2) := lambda_irSA2(r)
  end-if
end-do


! ============================================================
! CALCUL DES DEMANDES AGREGEES
! ============================================================

forall(r in R) do
  Lambda_r(r) := sum(i in I) lambda(r,i)
end-do

forall(i in I) do
  Lambda_i(i) := sum(r in R) lambda(r,i)
end-do


! ============================================================
! ! CALCUL DES Qstar PROPRES AU SCENARIO COURANT
! ============================================================

forall(d in D, r in R) do
  if Lambda_r(r) > 0 then
    Qstar(d,r) := sqrt(2.0 * T * maxlist(1,Ddr(r,d)) * Lambda_r(r) / Hr)
  else
    Qstar(d,r) := 0
  end-if
end-do


! ============================================================
! RECUPERATION DU Cdr FIXE CALCULE A PARTIR DE SA1 ET SA2
! ============================================================

if CapPercent=10 then
  CapIndex:=1
end-if

if CapPercent=40 then
  CapIndex:=2
end-if

if CapPercent=70 then
  CapIndex:=3
end-if

if CapPercent=100 then
  CapIndex:=4
end-if

Cdr := CdrFixed(MapID,CapIndex)

! ============================================================
! CALCUL DES Qdr CONTRAINTS
! ============================================================

forall(d in D, r in R) do
  if Lambda_r(r) > 0 then
    Qdr(d,r) := minlist(Cdr, Qstar(d,r))
  else
    Qdr(d,r) := Cdr
  end-if
end-do



nbContraints := 0
nbNonContraints := 0

forall(d in D, r in R) do
  if Qdr(d,r) = Cdr then
    nbContraints += 1
  else
    nbNonContraints += 1
  end-if
end-do

writeln("Nombre de couples contraints par la capacité = ", nbContraints)
writeln("Nombre de couples non contraints = ", nbNonContraints)
writeln("Pourcentage non contraint observé = ", 100.0*nbNonContraints/3600)


! ============================================================
! DÉCLARATION DES TYPES DE VARIABLES
! ============================================================

! Variables binaires
forall(d in D) do
  y(d) is_binary
end-do
forall(d in D, r in R) do
  x(d,r) is_binary
end-do

! Variables auxiliaires continues non-négatives
forall(d in D, i in I) do
  v1(d,i) >= 0
  v2(d,i) >= 0
end-do


 ! Nombre total d'expéditions traite par les DC
FrequencyMeasure :=
  sum(d in D, r in R)
    (Lambda_r(r) / Qdr(d,r)) * x(d,r)

! Composante fixe du cout d'ouverture
FixedOpenCost :=
  sum(d in D) Fd*y(d)

! Composante variable liee ? la fr?quence des expeditions
VariableOpenCost :=
  COpe * FrequencyMeasure

! Cout total d'ouverture
OpenCost := FixedOpenCost + VariableOpenCost
 	TransCost := sum(d in D, r in R)T * Ddr(r,d) * (Lambda_r(r)/Qdr(d,r)) * x(d,r)
 	StockRCost := sum(d in D, r in R, i in I) Hr * (Qdr(d,r)/2) * (lambda(r,i)/Lambda_r(r)) * x(d,r)
 	StockDCCost := sum(d in D, i in I)sqrt(2*Kd * hd) * v1(d,i)
 	SafetyRCost := sum(d in D, r in R, i in I)Hr * z_alpha * sigma(r,i) * sqrt(LT) * x(d,r)
 	SafetyDCCost := sum(d in D, i in I) hd * z_alpha * v2(d,i)


forall(d in D, i in I) do
  C_v1(d,i) := sum(r in R) lambda(r,i) * x(d,r)^2 <= v1(d,i)^2
end-do

forall(d in D, i in I) do
  C_v2(d,i) := sum(r in R) sigma(r,i)^2 * x(d,r)^2 <= v2(d,i)^2
end-do

forall(r in R) do
  C_assign(r) := sum(d in D) x(d,r) = 1
end-do

forall(d in D, r in R) do
  C_open(d,r) := x(d,r) <= y(d)
end-do

! ============================================================
! RÉSOLUTION ET AFFICHAGE
! ============================================================
writeln("Résolution en cours...")
minimize(OpenCost + TransCost + StockRCost + StockDCCost + SafetyRCost + SafetyDCCost)

TotalCostOut := getobjval
OpenCostOut := getsol(OpenCost)

FixedOpenCostOut := getsol(FixedOpenCost)
VariableOpenCostOut := getsol(VariableOpenCost)
FrequencyMeasureOut := getsol(FrequencyMeasure)

writeln("Valeur optimale : ", getobjval)
writeln("FixedOpenCost = ", FixedOpenCostOut)
writeln("VariableOpenCost = ", VariableOpenCostOut)
writeln("FrequencyMeasure = ", FrequencyMeasureOut)
writeln("COpe utilisé = ", COpe)
writeln("OpenCost = ", getsol(OpenCost))
writeln("TransCost = ", getsol(TransCost))
writeln("StockRCost = ", getsol(StockRCost))
writeln("StockDCCost = ", getsol(StockDCCost))
writeln("SafetyRCost = ", getsol(SafetyRCost))
writeln("SafetyDCCost = ", getsol(SafetyDCCost))


if OpenCostOut > 0 then
  writeln(
    "Part fixe observée = ",
    100.0 * FixedOpenCostOut / OpenCostOut,"%")
end-if

! ============================================================
! CALCUL KPI
! ============================================================

NbDC := 0

forall(d in D) do
  if getsol(y(d)) > 0.5 then
    NbDC += 1
  end-if
end-do

TotalLoad := 0
TotalCapacity := 0

forall(d in D,r in R) do
   if getsol(x(d,r)) > 0.5 then
      TotalLoad += Qdr(d,r)
      TotalCapacity += Cdr
   end-if
end-do

if TotalCapacity > 0 then
   LoadingRate := TotalLoad / TotalCapacity
else
   LoadingRate := 0
end-if

writeln("========================================")
writeln("RESEAU - SCENARIO ", SCENARIO)
writeln("========================================")

forall(d in D) do
  if getsol(y(d)) > 0.5 then
    write("DC ", d, " : ")

    forall(r in R) do
      if getsol(x(d,r)) > 0.5 then
        write(r, ", ")
      end-if
    end-do

    writeln
  end-if
end-do

TimeOut := getparam("XPRS_TIME")
GapOut := abs(getobjval - getparam("XPRS_BESTBOUND"))/abs(getobjval)


! ============================================================
! EXPORT EXCEL
! ============================================================
if not CalibrationMode then
	
if share = 80 then
  ShareIndex := 1
end-if

if share = 60 then
  ShareIndex := 2
end-if

if share = 40 then
  ShareIndex := 3
end-if


   if CapPercent=10 then CapIndex:=1
   	end-if
   if CapPercent=40 then CapIndex:=2 
   	end-if
   if CapPercent=70 then CapIndex:=3 
   	end-if
   if CapPercent=100 then CapIndex:=4 
   	end-if

if Fd=500 then
  FdIndex:=1
end-if

if Fd=1000 then
  FdIndex:=2
end-if

if Fd=2000 then
  FdIndex:=3
end-if

RowExport :=
4+(MapID-1)*108+(CapIndex-1)*27+(FdIndex-1)*9+(ShareIndex-1)*3+(SCENARIO-1)
  
! Sauvegarde résultats dans variables réelles
TotalCostOut := getobjval
OpenCostOut := getsol(OpenCost)

FixedOpenCostOut := getsol(FixedOpenCost)
VariableOpenCostOut := getsol(VariableOpenCost)
FrequencyMeasureOut := getsol(FrequencyMeasure)

TransCostOut := getsol(TransCost)
StockRCostOut := getsol(StockRCost)
StockDCCostOut := getsol(StockDCCost)
SafetyRCostOut := getsol(SafetyRCost)
SafetyDCCostOut := getsol(SafetyDCCost)

if OpenCostOut > 0 then
  ObservedFixedShareOut := FixedOpenCostOut / OpenCostOut
else
  ObservedFixedShareOut := 0
end-if

initializations to "mmsheet.excel:noindex;Resultats_Modele_Expeditions.xlsx"
	
MapID                 as ("[Resume$A" + strfmt(RowExport,0) + "]")
CapPercent            as ("[Resume$B" + strfmt(RowExport,0) + "]")
Fd                    as ("[Resume$C" + strfmt(RowExport,0) + "]")
TargetFixedShare      as ("[Resume$D" + strfmt(RowExport,0) + "]")
COpe                  as ("[Resume$E" + strfmt(RowExport,0) + "]")
SCENARIO              as ("[Resume$F" + strfmt(RowExport,0) + "]")

TotalCostOut          as ("[Resume$G" + strfmt(RowExport,0) + "]")
OpenCostOut           as ("[Resume$H" + strfmt(RowExport,0) + "]")
FixedOpenCostOut      as ("[Resume$I" + strfmt(RowExport,0) + "]")
VariableOpenCostOut   as ("[Resume$J" + strfmt(RowExport,0) + "]")
ObservedFixedShareOut as ("[Resume$K" + strfmt(RowExport,0) + "]")
FrequencyMeasureOut   as ("[Resume$L" + strfmt(RowExport,0) + "]")

TransCostOut          as ("[Resume$M" + strfmt(RowExport,0) + "]")
StockRCostOut         as ("[Resume$N" + strfmt(RowExport,0) + "]")
StockDCCostOut        as ("[Resume$O" + strfmt(RowExport,0) + "]")
SafetyRCostOut        as ("[Resume$P" + strfmt(RowExport,0) + "]")
SafetyDCCostOut       as ("[Resume$Q" + strfmt(RowExport,0) + "]")

NbDC                  as ("[Resume$R" + strfmt(RowExport,0) + "]")
LoadingRate           as ("[Resume$S" + strfmt(RowExport,0) + "]")
TimeOut               as ("[Resume$T" + strfmt(RowExport,0) + "]")
GapOut                as ("[Resume$U" + strfmt(RowExport,0) + "]")

end-initializations

writeln(
  "=> Resultats exportes dans Resultats_Modele_Expeditions.xlsx, ligne ",RowExport)

end-if

end-procedure


! ============================================================
! PRE-CALCUL DU Cdr FIXE POUR CHAQUE MAP ET CHAQUE Cap%
! A PARTIR DES Qstar DE SA1 ET SA2
! ============================================================

forall(m in 1..1) do

  MapID := m
  MapFile := "Maps" + strfmt(MapID,0) + "B_60R_50_250_500.xlsx"

  initializations from "mmsheet.excel:noindex;" + MapFile

    T as "[DonneesGenerales$N3]"
    Hr as "[DonneesGenerales$J2]"
    Ddr as "[Distances$B2:$BI61]"

    lambda_irSA1 as "[DonneesGenerales$C3:$C62]"
    lambda_irSA2 as "[DonneesGenerales$D3:$D62]"

  end-initializations

  forall(cp in {40}) do

    CapPercent := cp
    k := 0

    ! --------------------------------------------------------
    ! Qstar du scénario SA1 : 3600 valeurs
    ! --------------------------------------------------------

    forall(d in D, r in R) do

      Lambda_ref := lambda_irSA1(r)

      if Lambda_ref > 0 then

        k += 1

        QlistCap(k) :=
          sqrt(2.0 * T * maxlist(1,Ddr(r,d)) * Lambda_ref / Hr)
      end-if

    end-do

    ! --------------------------------------------------------
    ! Qstar du scénario SA2 : 3600 valeurs
    ! --------------------------------------------------------

    forall(d in D, r in R) do

      Lambda_ref := lambda_irSA2(r)

      if Lambda_ref > 0 then

        k += 1

        QlistCap(k) :=
          sqrt(2.0 * T * maxlist(1,Ddr(r,d)) * Lambda_ref / Hr)

      end-if

    end-do

    ! --------------------------------------------------------
    ! Tri des 7200 Qstar
    ! --------------------------------------------------------

    forall(ii in 1..7199) do
      forall(jj in ii+1..7200) do

        if QlistCap(ii) > QlistCap(jj) then

          temp := QlistCap(ii)
          QlistCap(ii) := QlistCap(jj)
          QlistCap(jj) := temp

        end-if

      end-do
    end-do

    ! --------------------------------------------------------
    ! Sélection du quantile correspondant : Cap%
    ! --------------------------------------------------------

    IndexCap := ceil(7200.0 * CapPercent / 100.0)

    if IndexCap < 1 then
      IndexCap := 1
    end-if

    if IndexCap > 7200 then
      IndexCap := 7200
    end-if

    if CapPercent=10 then
      CapIndex:=1
    end-if

    if CapPercent=40 then
      CapIndex:=2
    end-if

    if CapPercent=70 then
      CapIndex:=3
    end-if

    if CapPercent=100 then
      CapIndex:=4
    end-if

    CdrFixed(MapID,CapIndex) := QlistCap(IndexCap)

    writeln(
      "CdrFixed Map=",MapID,
      " Cap=",CapPercent," = ",CdrFixed(MapID,CapIndex))

  end-do
end-do


! ============================================================
! ANALYSE DE SENSIBILITE Fd - COpe
! TOUTES LES MAPS, Cap%, Fd ET PROPORTIONS
! ============================================================

forall(m in 1..1) do

  MapID := m

  forall(cp in {40}) do

    CapPercent := cp

    forall(fdv in {1000}) do

      Fd := fdv

      ! ======================================================
      ! CALIBRATION DE REFERENCE AVEC COpe = 0
      ! Cette calibration est faite une seule fois par
      ! combinaison Map - Cap% - Fd
      ! ======================================================

CalibrationMode := true
COpe := 0

      ! ----------------------------
      ! SA1 de référence
      ! ----------------------------

      SCENARIO := 1
      RunModel

BaselineFixedSA1 := FixedOpenCostOut
BaselineFrequencySA1 := FrequencyMeasureOut

      ! ----------------------------
      ! SA2 de référence
      ! ----------------------------

      SCENARIO := 2
      RunModel

BaselineFixedSA2 := FixedOpenCostOut
BaselineFrequencySA2 := FrequencyMeasureOut

      ! ----------------------------
      ! Agrégation SA1 + SA2
      ! ----------------------------

BaselineFixedTotal :=
  BaselineFixedSA1 + BaselineFixedSA2

BaselineFrequencyTotal :=
  BaselineFrequencySA1 + BaselineFrequencySA2

      ! ======================================================
      ! TEST DES DIFFERENTES PROPORTIONS
      ! ======================================================
      
      forall(sh in {60}) do
share := sh
TargetFixedShare := share / 100.0

writeln("share = ",share)
writeln("TargetFixedShare = ",TargetFixedShare)

        ! ----------------------------------------------------
        ! Calcul du coefficient COpe
        ! ----------------------------------------------------

if share = 100 then

  COpeCalculated := 0

elif share > 0 and BaselineFrequencyTotal > 0.000001 then

  COpeCalculated :=
    BaselineFixedTotal * (100.0 - share) / (share * BaselineFrequencyTotal)

else

  COpeCalculated := 0

end-if


        writeln("========================================")
        writeln("Map = ",MapID)
        writeln("Cap% = ",CapPercent)
        writeln("Fd = ",Fd)
        writeln(
          "Part fixe cible = ",
          100.0*TargetFixedShare,"%")
        writeln(
          "Coût fixe SA de référence = ",BaselineFixedTotal)
writeln(
  "Fréquence SA de référence = ",BaselineFrequencyTotal)
writeln("COpe calculé = ",COpeCalculated)
        writeln("========================================")

        ! ----------------------------------------------------
        ! Résolutions définitives
        ! M?me COpe dans SA1, SA2 et COOP
        ! ----------------------------------------------------

CalibrationMode := false
COpe := COpeCalculated

        forall(sc in 1..3) do

          SCENARIO := sc
          RunModel

        end-do

      end-do
    end-do
  end-do
end-do

end-model