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

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

  ! --- Paramètre ---
  T     : integer    
  Hr    : real      
  hd    : real    
  Kd    : integer     
  TotalLoad, TotalCapacity, LoadingRate  : 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 ---
  LT  : real    

  ! --- Lot sizes ---
  Qdr   : array(D,R) of real      
  Qstar : array(D,R) of real      ! permet de calculer Cap%
  QlistCap : array(1..7200) of real
  CdrFixed : array(1..5,1..4) of real
  Lambda_ref : 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, TransCost, StockRCost, StockDCCost, SafetyRCost, SafetyDCCost: linctr
  TotalCostOut, OpenCostOut, TransCostOut, StockRCostOut : real !export
  StockDCCostOut, SafetyRCostOut, SafetyDCCostOut : real
  
  
   ! --- Expériences ---
  MapID : integer
  NbDC : integer
  SCENARIO: integer
  
  
  ! --- Export Excel ---
  RowExport : integer
  CapIndex : 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
  C_OneLevel : array(D) of linctr
  C_CapacityK : array(D) of linctr
  
    !--- affichage GAP & TPS ---
  TimeOut, GapOut : real
  
    !--- Palier ---
  yK : array(D,K) of mpvar
  CapK : array(K) of real
  COpenK : array(K) of real
  NbP1, NbP2, NbP3, NbP4 : integer
  
 end-declarations
 
 
CapK(1) := 1000
CapK(2) := 2000
CapK(3) := 4000
CapK(4) := 8000

!COpenK(1) := 1000
!COpenK(2) := 1500
!COpenK(3) := 2250
!COpenK(4) := 3375

!---- Version 2 ----
COpenK(1) := 2000
COpenK(2) := 3000
COpenK(3) := 4500
COpenK(4) := 6750


procedure RunModel
	

   writeln("Run : Map=",MapID,
            " Cap=",CapPercent,
            " 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

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

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,kk in K) do
  yK(d,kk) 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

 
   ! --- Objectif ---
 	
 	OpenCost := sum(d in D,kk in K) COpenK(kk)*yK(d,kk)
 	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)


	! --- Contraintes ---
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) <= sum(kk in K) yK(d,kk)
end-do

! ----- Un seul palier par DC -----
forall(d in D) do
   C_OneLevel(d) :=
      sum(kk in K) yK(d,kk) <= 1
end-do

! ----- Capacité du DC -----
forall(d in D) do
   C_CapacityK(d) :=
      sum(r in R,i in I) lambda(r,i)*x(d,r) <= sum(kk in K) CapK(kk)*yK(d,kk)
end-do

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

writeln("Valeur optimale : ", getobjval)
writeln("OpenCost = ", getsol(OpenCost))
writeln("TransCost = ", getsol(TransCost))
writeln("StockRCost = ", getsol(StockRCost))
writeln("StockDCCost = ", getsol(StockDCCost))
writeln("SafetyRCost = ", getsol(SafetyRCost))
writeln("SafetyDCCost = ", getsol(SafetyDCCost))

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

NbDC := 0

forall(d in D) do
  if sum(kk in K) getsol(yK(d,kk)) > 0.5 then
    NbDC += 1
  end-if
end-do


! --- Capacité choisie ---
NbP1:=0
NbP2:=0
NbP3:=0
NbP4:=0

forall(d in D) do
 if getsol(yK(d,1))>0.5 then
    NbP1 +=1
 end-if

 if getsol(yK(d,2))>0.5 then
    NbP2 +=1
 end-if

 if getsol(yK(d,3))>0.5 then
    NbP3 +=1
 end-if
 if getsol(yK(d,4))>0.5 then
   NbP4 +=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


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


! ============================================================
! EXPORT EXCEL
! ============================================================

   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

 RowExport := 4 + (MapID-1)*12 + (CapIndex-1)*3 + (SCENARIO-1)

! Sauvegarde résultats dans variables réelles
TotalCostOut := getobjval
OpenCostOut := getsol(OpenCost)
TransCostOut := getsol(TransCost)
StockRCostOut := getsol(StockRCost)
StockDCCostOut := getsol(StockDCCost)
SafetyRCostOut := getsol(SafetyRCost)
SafetyDCCostOut := getsol(SafetyDCCost)

!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_1.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_2.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_3.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_4.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_5.xlsx"

!--- V2 ---
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_V2_1.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_V2_2.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_V2_3.xlsx"
!initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_V2_4.xlsx"
initializations to "mmsheet.excel:noindex;Resultats_Modele_Paliers_V2_5.xlsx"


MapID              as ("[Resume$A" + strfmt(RowExport,0) + "]")
CapPercent         as ("[Resume$B" + strfmt(RowExport,0) + "]")
SCENARIO           as ("[Resume$C" + strfmt(RowExport,0) + "]")

TotalCostOut       as ("[Resume$D" + strfmt(RowExport,0) + "]")
OpenCostOut        as ("[Resume$E" + strfmt(RowExport,0) + "]")
TransCostOut       as ("[Resume$F" + strfmt(RowExport,0) + "]")
StockRCostOut      as ("[Resume$G" + strfmt(RowExport,0) + "]")
StockDCCostOut     as ("[Resume$H" + strfmt(RowExport,0) + "]")
SafetyRCostOut     as ("[Resume$I" + strfmt(RowExport,0) + "]")
SafetyDCCostOut    as ("[Resume$J" + strfmt(RowExport,0) + "]")
NbDC               as ("[Resume$K" + strfmt(RowExport,0) + "]")

NbP1 as ("[Resume$L" + strfmt(RowExport,0) + "]")
NbP2 as ("[Resume$M" + strfmt(RowExport,0) + "]")
NbP3 as ("[Resume$N" + strfmt(RowExport,0) + "]")
NbP4 as ("[Resume$O" + strfmt(RowExport,0) + "]")
LoadingRate as ("[Resume$P"+strfmt(RowExport,0)+"]")

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

end-initializations

writeln("=> Resultats exportes dans Resultats_Modele_Paliers.xlsx")

writeln("===== RESEAU =====")

forall(d in D) do
  forall(kk in K) do
    if getsol(yK(d,kk)) > 0.5 then
      writeln("DC ", d, " - Palier ", kk, " :")
      forall(r in R) do
        if getsol(x(d,r)) > 0.5 then
          write(r, ", ")
        end-if
      end-do
      writeln
    end-if
  end-do
end-do

end-procedure

forall(m in 5..5) 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 {10, 40, 70, 100}) do
    CapPercent := cp

    k := 0

    ! Qstar pour SA1
    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 pour SA2
    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 valeurs
    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

    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


forall(m in 5..5) do
  forall(cp in {10, 40, 70, 100}) do
    forall(sc in 1..3) do

      MapID := m
      CapPercent := cp
      SCENARIO := sc

      RunModel
    end-do
  end-do
end-do

end-model