Decision Optimization

Decision Optimization

Delivers prescriptive analytics capabilities and decision intelligence to improve decision-making.

 View Only

  • 1.  User functions in model. Roster model

    Posted Mon June 21, 2021 01:50 PM
    Hi,

    I've been struggling for some time with this model. It's a rostering problem. I have a working model in minizinc (using cplex as solver) but I am not satisfied with the performance so I have decided to implement it directly in cplex to see if there's an improvement. The objective of the model is to calculate a valid roster (CalculatedRoster) meeting some criteria.

    array[Staff, 1..31] of var TypeOfShift: CalculatedRoster;

    The model is as followed:

    - The staff could work different type of shifts (free (l), early-morning(m1), morning(m), evening(t), night(n), oncall-morning(im), oncall-evening(it), oncall-night(ino), office(o)) each of them defined with a starting time and duration arrays (startShift, durationShift). In minizinc you can use the enum type that makes the code quite readable but it should not be a problem to assign to each shift a number (0 to 9) for cplex. All time are in minutes.

    enum TypeOfShift = {l,m1,m,t,n,im,it,ino,o}
    array[TypeOfShift ] of int: startShift=[600, 405, 450, 885, 1320, 390, 825, 1260, 600];
    array[TypeOfShift ] of int: durationShift=[0, 435, 435, 435, 570, 495, 495, 630, 360];

    The free shift (l) is defined with a starting time of 600 and duration 0 so to meet the 12h resting period between shifts that I will explain later.
    An example of a roster for 3 people and 31 days:
    l, l, l, t, t, t, it, t, n, l, l, l, l, l, l, l, l, l, l, l, im, m1, t, t, ino, l, l, l, m1, m, t
    n, l, l, l, l, l, l, l, o, o, l, l, l, l, l, l, l, l, l, l, it, t, t, ino, n, l, l, l, t, t, t
    m, l, l, l, m, m, im, im, m, l, t, t, t, t, n, l, l, l, l, l, m, m, m, n, n, l, l, l, l, l, l

    - There are different criteria that need to be complied for a roster to be valid: one of them is that there has to be a 12h (720 minutes) between the end of a shift and the start of the following. This constraint should be easy to implement in cplex, in minizinc is like this:

    constraint forall(i in Staff, d in 1..30)
    ((CalculatedRoster[i,d+1] != l ) -> (1440*(d-1)+startShift[CalculatedRoster[i,d]] + durationShift[CalculatedRoster[i,d]] + 720 <= 1440*d+startShift[CalculatedRoster[i,d+1]]));
    - There are special resting periods after some shifts:
      1. After one night (48h)
    2. After two consecutive nights (54h = 3240 minutes). nn or nln are both considered consecutive nights.
    3. After a rotation of 5 shifts (48h = 2880 minutes) not including on those shifts the office (o)
    4. After a rotation of 6 shifts (54h) not including on those shifts the office (o). 6 is as well the maximum number of shifts in a rotation (not including office)
    5. Maximum of 2 oncall shifts per rotation
    1 and 2 should be easy to implement in cplex as it's a straight forward syntax translation from minizinc. The problem is with 3 and 4. The way I did it in minizinc was with a state machine definition (it's called regular constraint), however I think this is not available in cplex. I implemented as well with the below code which uses a user function (RestTime) (not available either in cpex). Basically, the idea is to detect the first day where a resting period of >=54h (2) or  >= 48h (1) is available and set the rest of the days to 0. rest = 1440*(nextShift-day) + startShift[CalculatedRoster[i,nextShift]] - (startShift[CalculatedRoster[i,dia]] + durationShift[CalculatedRoster[i,day]] calculates the difference in minutes between the end of a shift and the start of the following.

    %Rotation
    array[Staff, 0..31] of var 0..2: Rotation = array2d(Staff, 0..31, [if(d == 31) \/ d == 0 then 2 elseif CalculatedRoster[i,d] == l then 0 else RestTime(i,d) endif | i in Staff, d in 0..31]);
    %RestTime
    function var 0..2: RestTime(Staff: i, int: day) =
    let {var day+1..31: nextShift = min([d | d in day+1..31 where CalculatedRoster[i,d] != l]++[31]),
    var int: rest = 1440*(nextShift-day) + startShift[CalculatedRoster[i,nextShift]] - (startShift[CalculatedRoster[i,dia]] + durationShift[CalculatedRoster[i,day]]) } in
    (if rest >= 3240 then 2 elseif (rest >= 2880) then 1 else 0 endif);

    For the example above the rotation array would be:
    2;0;0;0;0;0;0;0;0;2;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;2;0;0;0;0;0;2;
    2;2;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;2;0;0;0;0;0;2;
    2;2;0;0;0;0;0;0;0;1;0;0;0;0;0;2;0;0;0;0;0;0;0;0;0;2;0;0;0;0;0;2;

    Once I have this array, I can work the 3, 4 and 5 criteria (and more). The basic idea is to sum shifts between Rotation that are not 0.

    constraint forall (i in Staff, d in 0..29 where Rotation[i,d] != 0)
    %d is the beginning of the rotation and d1 the end of the rotation
    (let {var int: d1 = min([dx | dx in d+1..30 where Rotation[i,dx] != 0 ]++[31])} in
    %maximum number of shifts 6
    (sum([1 | dx2 in d+1..d1 where CalculatedRoster[i,dx2] != l /\ CalculatedRoster[i,dx2] != o]) <= 6) /\
    %maximum number of morning shifts in a rotation is 5
    (sum([1 | dx2 in d+1..d1 where CalculatedRoster[i,dx2] == m \/ CalculatedRoster[i,dx2] == m1]) <= 5) /\
    %54h rest after 6 days (the 48h rest after 5 days is meet by default with previous and next constrain)
    (sum([1 | dx2 in d+1..d1 where CalculatedRoster[i,dx2] != l /\ TurnCalculatedRosterroCalculado[i,dx2] != o]) == 6 -> Rotation[i,d1] == 2) /\
    %maximum number of oncall shifts in a rotation is 2
    (sum([1 | dx2 in d+1..d1 where CalculatedRoster[i,dx2] == im \/ CalculatedRoster[i,dx2] == it \/ CalculatedRoster[i,dx2] == ino]) <=2));
    again, I am using the minizinc possibility of defining a local variable inside a scope (let {var int: d1 = ...... ) that is not available in cplex.

    Sorry for the long post but I have tried to explain the complexity of the model good enough, even though I have simplified some parts. I have read other roster models but the complexity of the criteria is different. I know cplex has special functions for schedulling but I have not been able to fit them in my model. May be, there's another way to aproach this model. Any help would be much appreciated

    Cheers

    ------------------------------
    javier rodrigo
    ------------------------------

    #DecisionOptimization


  • 2.  RE: User functions in model. Roster model

    Posted Sat June 26, 2021 01:45 PM
    Edited by System Admin Fri January 20, 2023 04:19 PM
    I manage to implement the whole model in CP but when I activate the rotations constraints (counting the shifts in between rests), the engine never finish the extraction. I believe there are just too many constraints and variables (70000 and 9000 respectively). I will try to focus my problem. This is in short what I have now:

    tuple ShiftPattern{
    key string Type;
    int Equiv; //integer equivalent of the shift Type
    int initTime; //starting time of the shift in minutes
    int Duration; //duration of the shift in minutes
    }
    {ShiftPattern} TypeOfShift={ 
    <"l",0,0,0>, <"m1",1,405,435>, <"m",2,450,435>, <"im",3,390,495>, <"o",4,600,360>, <"t",5,885,435>, <"it",6,825,495>,
    <"n",7,1320,570>, <"ino",8,1260,630>};

    dvar int CalculatedRoster[Staff][1..31][TypeOfShift] in 0..1

    As I said there are several criterias to comply but the rotation ones are the hardest. To calculate which days are in each rotation in a month period I need to calculate the resting periods >48h and >54h. I do it like this (i know, not easy):

    dexpr int rest[i in Staff][d in 1..31] = (CalculatedRoster[i][d][<"l">] == 0)*(100000 - maxl(0,max(dia in d+1..31, t in TypeOfShift: t.Tipo != "l")(100000-(1440*(dia-d)+t.initTime- (sum(t2 in TypeOfShift: t2.Tipo != "l") (CalculatedRoster[i][d][t2]*(t2.initTime+ t2.Duration)))))*(CalculatedRoster[i][dia][t]==1)));

    for the following example
    l, l, l, t, t, t, it, t, n, l, l, l, l, l, l, l, l, l, l, l, im, m1, t, t, ino, l, l, l, m1, m, t
    n, l, l, l, l, l, l, l, o, o, l, l, l, l, l, l, l, l, l, l, it, t, t, ino, n, l, l, l, t, t, t
    m, l, l, l, m, m, im, im, m, l, t, t, t, t, n, l, l, l, l, l, m, m, m, n, n, l, l, l, l, l, l

    rest is:
    [0 0 0 1005 1005 945 1005 1440 15780 0 0 0 0 0 0 0 0 0 0 0 960 1485 1005 1380 4275 0 0 0 1050 1440 100000]
    [8790 0 0 0 0 0 0 1080 17145 0 0 0 0 0 0 0 0 0 0 0 1005 1005 1380 870 4755 0 0 0 1005 1005 100000]
    [5325 0 0 0 1005 945 945 1005 2880 0 1005 1005 1005 1440 7200 0 0 0 0 0 1005 1005 1875 870 100000 0 0 0 0 0 0]

    and then using 
    dexpr int Rotations[i in Staff][d in 0..31] = (d == 31) || (d == 0) ? 2 : rest[i][d] >= 3240 ? 2 : rest[i][d] >= 2880 ? 1 : 0;
    I get this:
    [2 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2]
    [2 2 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2]
    [2 2 0 0 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2]
    where 2 means >=54h rest, 1 means >=48hs and 0 means <48hours

    As these are variables, it has to be computed many times, so it has to be fast. Now what I need based on the Rotations and the CalculatedRoster is to make the constraints of number of shifts in each rotation. This is where I am currently stuck. I need to count the number of shifts different than "l" or "o" from CalculatedRoster in between each 1/2 of Rotations . This number has to be <= 5 is the number is 1 and <=6 if the number is 2. Any ideas?

    When I run a very simplified model without the rotation constraints I get almost inmediately that there's no solution (which is correct in this specific test). However, then it starts trying until reaches the same conclussion... no solution. Why is that? (see image)



    Any idea if trying to implement the model as schedulling problem, I mean using interval variables and all these features, would make a different? I am not familiar at all with that so it would take me a great effort to implement it...

    I have seen the nurse shift example and that's why I decided to implemented that way but the model is much more bigger (31 workers, 31 days and 9 type of shifts per day)

    thanks

    ------------------------------
    javier rodrigo
    ------------------------------

    Original Message


  • 3.  RE: User functions in model. Roster model

    Posted Mon June 28, 2021 02:49 AM
    Hi,

    1) Can you post your entire .mod and .dat so that other users could try ?

    2)

    Instead of

    dexpr int rest[i in Staff][d in 1..31] = (CalculatedRoster[i][d][<"l">] == 0)*(100000 - maxl(0,max(dia in d+1..31, t in TypeOfShift: t.Tipo != "l")(100000-(1440*(dia-d)+t.initTime- (sum(t2 in TypeOfShift: t2.Tipo != "l") (CalculatedRoster[i][d][t2]*(t2.initTime+ t2.Duration)))))*(CalculatedRoster[i][dia][t]==1)));​


    could you try to write that as a decision variable dvar ?

    regards



    ------------------------------
    [Alex] [Fleischer]
    [EMEA CPLEX Optimization Technical Sales]
    [IBM]
    ------------------------------

    Original Message


  • 4.  RE: User functions in model. Roster model

    Posted Tue June 29, 2021 10:55 AM
    but rest is a direct calculation from CalculatedRoster. Is it faster for the solver to find a solution for rest (as a dvar) rather than to calculate it? Nevertheless, I will try to put it as  a dvar

    What about expressing the model with interval variables? Is it worthy?

    ------------------------------
    javier rodrigo
    ------------------------------

    Original Message


  • 5.  RE: User functions in model. Roster model

    Posted Tue July 06, 2021 08:53 AM
    Hi,

    This is the full model (sorry for the spanish naming). Even that I am using sol.setValue to give an initial value to the solution and that the solution and the initial values are always very close (normally only a few shifts change between then). It is just taking very long. I have even reduce the total amount of days from 31 to a smaller window, but it is still taking long (>5minutes). I have tried with different search methods but it was just a trial and error as I don´t understand what they do. Any ideas on how to improve the performance?

    regards

    ------------------------------
    javier rodrigo
    ------------------------------

    Original Message


  • 6.  RE: User functions in model. Roster model

    Posted Tue July 06, 2021 11:00 AM
    Hi

    I would rather rely on MIP and write the following which gives

    Multi-objective solve log . . .

Index  Priority  Blend          Objective      Nodes  Time (sec.)  DetTime (ticks)
    1         3      1   2,0000000000e+00          0       131,38        164980,36
    2         2      1   0,0000000000e+00          0       135,06       1641846,59
    3         1      1   4,0000000000e+00          0       140,41       1642106,66
    4         0      1   0,0000000000e+00          0       154,25       1648991,16​


    in the engine log and 



    DR l 8 0 1
DR n 8 1 0
RC l 10 0 1
RC m 10 1 0
MW l l t t t t it n l l l l l l l l l l l l im m1 t t ino l l l m1 m t 
RS n l l l l l l l l l l l l l l l l l l l it t t ino n l l l t t t 
AL m l l m m m im im l l t t t t n l l l l l m m m n n l l l l l l 
JP l l m l m m m ino n l l l m m1 t ino n l l l m t t t n l l l m m1 t 
SG l l n l l t t t ino n l l l it t n n l l l l m it t t n l l l t it 
CE l n l l l im im t t n l l l im m t t n l l l l l l l l l l l m m 
JG t n l l l m m1 m m1 o o l l l m t t n l l l m m t t n l l l im im 
DV t n n l l it t t t n l l l t t t ino n l l l im im m t n l l l t t 
AY m l l l ino n l l it ino n l l l m t t t n l l l m m m1 ino n l l l m 
SL l l l l l l l l l l l l l l im m1 t ino n l l l m1 m m t n l l l m1 
RJ l m1 n l l t l l l l l l l l o o o it n l l l o o o m1 n l l l im 
HL m l n l l l m m t it n l l l m m t t n l l l m m it t n l l l m 
SH t t l n l l l m m t n n l l l l l l l l l l l it t t ino n l l l 
SM l t n n l l l im m m ino n l l l m m t t n l l l o o it t n l l l 
EG m1 t l l o o l l o o o o o l l m im t t n l l o o o l o l l o o 
VC l l l l l l l it t t n n l l l m m1 t ino n l l l m1 t t t n l l l 
AM m m t l n l l l m1 m1 t t n l l l m m m t n l l l m m it ino n l l 
JK l l l l l l l l o m m1 m n l l l m m1 im im n l l l l l l l l l l 
JL t t t n n l l l t t t l ino l l l it t t ino n l l l m t t t n l l 
CH l m m1 t n n l l l o o t t n l l l m it t t n l l l m t t ino n l 
LB l m t t n n l l l m it t t n l l l m m1 it ino n l l l l l l l l l 
EE n l l l l l l l l t t t n n l l l im m m im n l l l m m1 m1 n n l 
CF l l m m t ino n l l l m m m m n l l l l l l l l l l l m it t t ino 
MC l l l l l l l l l l t t t t ino l l l t t it t n l l l t t t ino n 
DR l t t t t t l n l l m it t t n l l l m m1 t ino n l l l l l l l l 
RC l l m t t t n l l m l l l l l l l l m it t t ino l l l m t t t n 
GG l l m m m n n l l l m t it ino n l l l m t t it n l l l m m t n n 
JJ l l l m m t t n l l l im m t t n l l l m m m o n l l l m im it t 
BB l l l m1 m1 t it n l l l l l l l l l l l m m t t n l l l m m t it 
SC l l o o o m1 ino n l l l m1 im t it n l l l m m1 t it n l l l l l l l 
PT l l l m it t it n l l l m m1 m m1 n l l l l l l l l l l l m m m m 
​

    in the scripting log after 10 minutes
    /* Modelo para buscar qué cta puede realizar un determinado turno bien por asignación directa o mediante cambios. Con matrices de transición */
/*
Si el modelo se utiliza para buscar un cambio (BusquedaDeUnCambio=1), CtaActual contendrá el valor del cta que busca el cambio (CtaActual=X)
Si el modelo se utiliza para buscar hhee (BusquedaDeUnCambio=0), se puede llamar con CtaActual=X siendo X el cta al que le quiero asignar las hhee. El .bat es el que llama al modelo tantas veces como ctas hay (aptos según excel)
dvar int TurneroCalculado[Plantilla][1..NumDiasTurnero][turnos] in 0..1; MAS RAPIDO AL EJECUTAR LAS CONSTRAINT DE NN, NLN QUE LA V2
definir Ciclos como dvar en lugar de dexpr. Funciona más rápido que minizinc (se podrían implementar en minizinc algunas de estas mejoras)
Mas lento que el minizin con matrices
*/
//using CP;



string CtaActual="TODOS"; //SI SE LLAMA AL MODELO DESDE EXCEL QUITAR LA ASIGNACION DEL VALOR. SI SE LLAMA DESDE CPLEX PONER la inicial del cta. "TODOS" equivale a CtaActual=0

//Definción Parámetros de entrada
int NumCtas=...; //Número de ctas que se pasa en la matriz del turnero
{string} Plantilla = ...;  //Listado de iniciales de los ctas
int NumDiasTurnero=...; //Número de días que se pasan en la matriz de turnero, inamovible y duración
int DiasTurnero[1..NumDiasTurnero]=...; //Días a los que corresponde el turnero
string Turnero[Plantilla, 1..NumDiasTurnero]=...;  //Turnero de trabajo
int Inamovibles[Plantilla, 1..NumDiasTurnero]=...;  //Matriz de turnos que se pueden (0) y no se puede mover (1)
int FechaTurnoABuscar[1..5]=...;  // Días que se quiere buscar de hhee (índice del array de DiasTurnero) o cambios. Si es 0 implica que no se busca nada
string TurnoABuscar[1..5]=...;   // Turno que se quiere buscar de hhee o cambios
int DiaActual=...;    // Día actual para impedir que trate de hacer cambios del pasado (índice del array de DiasTurnero)
int BusquedaDeUnCambio=...;  // Estoy buscando un cambio (no cubrir unas hhee)
int DesprogramarNoche3mas1=...; // Desprogramar noches de 3+1 solo para buscar hhee
int DesprogramarImag=...; // Desprogramar imaginarias

//busca hasta FranjaSupDiasTurn días por encima de la fecha de cambio/hhee. Para limitar el espacion de búsqueda y acelerar
int FranjaSupDiasTurn=10;
int DiaLimite = max(d in 1..5) FechaTurnoABuscar[d] + FranjaSupDiasTurn;

//Definiciones de los tipos de turnos: m, im... (l a lo mejor no hace falta)
tuple PatronTurno {
  key string Tipo;
  int Equiv; //equivalente numérico para trabajar en el modelo
  int HoraInicio; //horarios de inicio de los distintos tipos de turnos (en minutos)
  int Duracion; //duración de los distintos tipos de turnos. Las RJs y Oficinas no es lo exacto (se podría poner la duración exacta defininiendo todos los posibles turnos de Oficina y RJs) (en minutos)
}
{PatronTurno} tipoDeTurno={ //Los tipos de turnos que admite el modelo (inicial) (**************************************PROBAR TAMBIEN A DEFINIR LOS TIEMPOS EN HORAS COMO FLOAT)
<"l",0,0,0>, <"m1",1,405,435>, <"m",2,450,435>, <"im",3,390,495>, <"o",4,600,200>, <"t",5,885,435>, <"it",6,825,495>, 
<"n",7,1320,570>, <"ino",8,1260,630>};

tuple Turno {
  key string Tipo;
  int Equiv; //equivalente numérico para trabajar en el modelo
  key int Dia; //día del turnero donde se realiza el turno
  int AbsHoraInicio; //horarios de inicio de los distintos tipos de turnos (en minutos)
  int AbsHoraFin; //duración de los distintos tipos de turnos. Las RJs y Oficinas no es lo exacto (se podría poner la duración exacta defininiendo todos los posibles turnos de Oficina y RJs) (en minutos)
}
//Crea todos los turnos del modelo
{Turno} turnos = {<t.Tipo, t.Equiv, d, t.HoraInicio+1440*(d-1), t.HoraInicio+1440*(d-1)+t.Duracion> | d in 1..NumDiasTurnero, t in tipoDeTurno};

//{string} TipoDeTurnoStr = {"l", "m1", "m", "t", "n", "im", "it", "ino", "o"};  //Los tipos de turnos que admite el modelo
//int HoraInicioTipoDeTurno[0..card(TipoDeTurno)-1] = [600, 405, 450, 390, 600, 885,  825, 1320, 1260]; //horarios de inicio de los distintos tipos de turnos (en minutos) indexado con el valor del turno
//int DuracionTipoDeTurno[0..card(TipoDeTurno)-1] = [0, 435, 435, 495, 360, 435, 495, 570, 630]; //duración de los distintos tipos de turnos. indexado con el valor del turno

//TuneroAct y TuneroActualizado (pasando de letras a turnos para trabajar en el modelo)
//06/06/21 Actualizo turnero (creo TurneroActualizado): con todas las imaginarias anteriores al DiaActual las convierto en días libres (si una imaginaria fue activada debería venir como un turno (m,t,n) de Excel). 
//Al abrir el resultado en excel se reconvierten los "l" en sus imaginarias
string TurneroAct[i in Plantilla][d in 1..NumDiasTurnero] = 
	((Turnero[i][d] == "im" || Turnero[i][d] == "it" || Turnero[i][d] == "ino") && (d < DiaActual)) ? "l" : Turnero[i][d];
int TurneroActualizado[i in Plantilla][d in 1..NumDiasTurnero][t in tipoDeTurno] = TurneroAct[i][d] == t.Tipo ? 1 : 0 ; //contiene 1 en el turno que se hace (libre no está definido, sería todo a cero)

int InamoviblesNueva[i in Plantilla][d in 1..NumDiasTurnero]; //define una nueva matriz inamovible
execute{
	//cp.param.timelimit=TiempoHastaEncontrarUnaSolucion;  //fija un límite temporal
	//cp.param.LogVerbosity= "Quiet";	//evita que vaya escribiendo en el log file las iteraciones para encontrar soluciones	
	//06/06/21 Matrices inamovibles: Matriz de turnos que se pueden (0) y no se puede mover (1) para fijar que ctas pueden modificar su turno
	//Crea la matriz InamoviblesNueva igual que Inamovibles. Además si buscamos un CtaActual != TODOS añade 1 en la fecha correspondiente para que se le asigne y no se mueva el turno que buscamos. Además si es d < DiaActual lo pone todo a 1 
	for (var i in Plantilla){
		for(var d=1; d<=NumDiasTurnero;d++){
		  	if (CtaActual == "TODOS"){	
				InamoviblesNueva[i][d] = Inamovibles[i][d]
 			}
 			else {
 			  	if ((i == CtaActual) && ((d==FechaTurnoABuscar[1]) || (d==FechaTurnoABuscar[2]) || (d==FechaTurnoABuscar[3]) || (d==FechaTurnoABuscar[4]) || (d==FechaTurnoABuscar[5])))
 			  		InamoviblesNueva[i][d]=1
 			  	else {
 			  		InamoviblesNueva[i][d] = Inamovibles[i][d]
     			} 			  	
 			}			
		}	
	}
}

// Calculo de CRs del turnero
int CRs[d in 1..NumDiasTurnero][t in tipoDeTurno] = sum(i in Plantilla) (TurneroActualizado[i][d][t]); // CRs para las m,t,n, im, it, ino, oficinas... En el índice 0 van las l, en el 1 van las m1
int CRsNuevas[1..NumDiasTurnero][tipoDeTurno]; //Si se busca cambio se mantienen las CRs originales; si no son las CRs del turnero original y añadido el turno a buscar
//int CRsNuevas[1..NumDiasTurnero][tipoDeTurno] = CRs; //USAR ESTA DEFINICION PARA chequeo CRsNuevas=CRs

execute{ // descomentar esta si estamos haciendo uso normal, es decir, si no es para debugging
	for (var d=1; d<=NumDiasTurnero; d++){
		for(var t in tipoDeTurno){
			if (BusquedaDeUnCambio == 0){	//buscamos hhee
				if (((d==FechaTurnoABuscar[1]) && (t.Tipo==TurnoABuscar[1])) || ((d==FechaTurnoABuscar[2]) && (t.Tipo==TurnoABuscar[2])) 
               		|| ((d==FechaTurnoABuscar[3]) && (t.Tipo==TurnoABuscar[3]))  || ((d==FechaTurnoABuscar[4]) && (t.Tipo==TurnoABuscar[4]))
               		|| ((d==FechaTurnoABuscar[5]) && (t.Tipo==TurnoABuscar[5]))) {
					CRsNuevas[d][t] = CRs[d][t]+1
  				}				
				else CRsNuevas[d][t] = CRs[d][t]
 			}
 			else CRsNuevas[d][t] = CRs[d][t]		  			
		}
	}			
}

//Turnos seguidos incompatibles
//["t", "m"], ["it", "m"], ["t", "m1"], ["it", "m1"], ["t", "im"], ["it", "im"], ["t", "o"], 
//			["it", "o" ], ["n", "m"], ["ino", "m"], ["n", "m1"], ["ino", "m1"], ["n", "im"], ["ino", "im"], 
//			["n", "o"], ["ino", "o"], ["n", "t"], ["ino", "t"], ["n", "it"], ["ino", "it"]];
//Turnos consecutivos incompatibles: los que no respeta la regla de las 12h entre ellos 
tuple turnoIncompatible{
  string t_i;
  string t_f;
}
//PLANTEARLO COMO ARRAY POR SI SE PUDIERA ACELERAR ESA CONSTRAING ************************************************************************************************
{turnoIncompatible} turnosIncompatibles = {<t1.Tipo, t2.Tipo> |  t1,t2 in turnos: t1.Tipo != "l" && t2.Tipo != "l" && (t1.AbsHoraFin + 720 > t2.AbsHoraInicio) && t1.Dia == 1 && t2.Dia == 2};

//Combinación de turnos que consiguen descansos de 48h o más															
tuple rompeCiclo3D{
  string t1;
  string t2;
  string t3;
  int desc;
}
{rompeCiclo3D} t_rompeCiclo3D = {<t1.Tipo, t2.Tipo, t3.Tipo, ((t3.AbsHoraInicio-t1.AbsHoraFin) >= 3240 || t3.Tipo == "l" ? 2 : 1) > 
								| t1,t2,t3 in turnos: t1.Dia == 1 && t2.Dia == 2 && t3.Dia == 3 && t1.Tipo != "l" && t1.Tipo != "n" && t1.Tipo != "ino" && t2.Tipo == "l" && //t1 no puede ser n o in porque necesitamos un cuarto día para poder ver si hay >48h descanso
							((t1.AbsHoraFin + 2880 <= t3.AbsHoraInicio) || t3.Tipo == "l") };
							
tuple rompeCiclo4D{
  string t1;
  string t2;
  string t3;
  string t4;
  int desc;
}
{rompeCiclo4D} t_rompeCiclo4D = {<t1.Tipo, t2.Tipo, t3.Tipo, t4.Tipo, ((t4.AbsHoraInicio-t1.AbsHoraFin) >= 3240 || t4.Tipo == "l" ? 2 : 1)> | t1,t2,t3,t4 in turnos: 
				t1.Dia == 1 && t2.Dia == 2 && t3.Dia == 3 && t4.Dia == 4 && t1.Tipo != "l" && t1.Equiv >= 7 && t2.Tipo == "l" && t3.Tipo == "l" && //solo miro los casos de dias con t1 sea n o ino
							((t1.AbsHoraFin + 2880 <= t4.AbsHoraInicio) || t4.Tipo == "l") };
							

// Variables OJO CHEQUEAR SI ESTAN BIEN DEFINIDAS ****************************************
// CONSIDERAR CREAR UN DVAR TURNEROCALCULADO COMO EN MINIZINC (PONIENDO INT EN LUGAR DE TURNOS) PARA ACELERARAR ALGUNAS RESTRICCIONES
dvar int TurneroCalculado[Plantilla][1..NumDiasTurnero][tipoDeTurno] in 0..1;
//dexpr int TurneroCalculado2[i in Plantilla][t in turnos] = TurneroCalculado[i][t.Dia][<t.Tipo>] ;
//dexpr int TurneroCalculadoDias[i in Plantilla][d in 1..NumDiasTurnero] = max(t in turnos : t.Dia == d) TurneroCalculado[i][t]*t.Equiv;
dexpr int Cambios = sum(i in Plantilla, d in 1..NumDiasTurnero, t in tipoDeTurno) (TurneroActualizado[i][d][t] != TurneroCalculado[i][d][t]); // número de cambios que se ha hecho con respecto al turnero original.
dexpr int DesprogOficina = sum(i in Plantilla, d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "o") (TurneroActualizado[i][d][t] != TurneroCalculado[i][d][t]); // desprogramación de oficina
dexpr int DesprogImag = sum(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "im" || t.Tipo == "it" || t.Tipo == "ino") CRsNuevas[d][t] - 
						sum(i in Plantilla, d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "im" || t.Tipo == "it" || t.Tipo == "ino") (TurneroCalculado[i][d][t]); // número de desprogramación de imaginarias o no activación
dexpr int DesprogNoche = sum(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "n") CRsNuevas[d][t] - sum(i in Plantilla,  d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "n") (TurneroCalculado[i][d][t]); // //Desprogramaciones de noches de 3+1
dexpr int NumCtasAfectados =  sum (i in Plantilla) ((sum(d in 1..NumDiasTurnero, t in tipoDeTurno) (TurneroActualizado[i][d][t] != TurneroCalculado[i][d][t])) >= 1); // Número de ctas que hay que mover para colocar el nuevo turno						

//Calculo de donde caen los ciclos
//DICEN QUE ES MAS RAPIDO CONSTRAINTS QUE DEXPR
//tiempo inicio del siguiente día que tiene un turno != libre. Si no existe coge por defecto NumDiasTurnero  
//dexpr int Descanso[i in Plantilla][d in 1..NumDiasTurnero] = (TurneroCalculado[i][d][<"l">] == 0)*(100000 - maxl(0,max(dia in d+1..NumDiasTurnero, t in tipoDeTurno : t.Tipo != "l") 
//											(100000-(1440*(dia-d)+t.HoraInicio - (sum(t2 in tipoDeTurno : t2.Tipo != "l") (TurneroCalculado[i][d][t2]*(t2.HoraInicio + t2.Duracion)))))*(TurneroCalculado[i][dia][t]==1)));
//Ciclos: array que contiene los descansos de 48h o más
//vale 2 54h, 1 48h, 0 < 48h. El número está en el último turno del ciclo anterior y en el primer día del siguiente ciclo. Dia 0 y NumDiasTurnero es 2

//dexpr int Ciclos[i in Plantilla][d in 0..NumDiasTurnero] = (d == NumDiasTurnero) || (d == 0) ? 2 : Descanso[i][d] >= 3240 ? 2 : Descanso[i][d] >= 2880 ? 1 : 0;
dvar int Ciclos [i in Plantilla][d in 0..NumDiasTurnero] in 0..2;


//minimize NumCtasAfectados;
minimize staticLex(NumCtasAfectados, DesprogOficina, Cambios+DesprogImag, DesprogNoche);
subject to {
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//Para poder chequear si los constraints funcionan sobre el turnero que metemos. Además marcar como comentario la definción de CRsNuevas de arriba y poner CtaActual="TODOS". 
// Así chequea si el turnero que se ha metido cumple 1001  
// QUITAR SI NO ES PARA EL CHEQUEO
//forall(i in Plantilla, d in 1..NumDiasTurnero, t in tipoDeTurno)
//	(TurneroCalculado[i][d][t]==TurneroActualizado[i][d][t]);

//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

	//Si CtaActual=TODOS busca la mejor solución entre todos los ctas. Si CtaActual!=TODOS busca los cambios para que ese cta haga el turno
	if (CtaActual != "TODOS"){
		forall(f in 1..5 : FechaTurnoABuscar[f] != 0, t in tipoDeTurno : t.Tipo == TurnoABuscar[f]) 
			(TurneroCalculado[CtaActual][FechaTurnoABuscar[f]][t] == 1);	  
	} 
	 
	//Un único turno por cta
	forall(i in Plantilla, d in 1..NumDiasTurnero)
		((sum(t in tipoDeTurno) TurneroCalculado[i][d][t]) == 1);
		
              
	//Lo que ha pasado no se puede mover
	forall(i in Plantilla, d in 1..NumDiasTurnero : d<=DiaActual || d>DiaLimite, t in tipoDeTurno)
    	(TurneroCalculado[i][d][t]==TurneroActualizado[i][d][t]);
  
 
 	//Lo que es inamovible no se mueve (tiene en cuenta cuando se fuerza a que un turno sea inamovible para que lo haga un determinado cta (en este caso InamoviblesNueva[i,d] será 1 y será distinto de Inamovibles[i,d])
	forall(i in Plantilla, d in 1..NumDiasTurnero : InamoviblesNueva[i][d]==1 && InamoviblesNueva[i][d]==Inamovibles[i][d], t in tipoDeTurno)
		(TurneroCalculado[i][d][t]==TurneroActualizado[i][d][t]);
  
	//Las v y oficinas no se puede mover ni crear (la condición para las v está dentro de la de CRs e inamovibles)
	forall(i in Plantilla, d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "o")
		(TurneroActualizado[i][d][t] == 0) => (TurneroCalculado[i][d][t] == 0);
	
	//Las CRs hay que mantenerlas obligatoriamente para m, m1 y t. Si fuerzo también a manterner las CRs de v como son inamovibles entonces estoy forzando a que no se generen mas v
	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "m")
		(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);
		
	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "m1")
		(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);
		
	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "t")
		(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);
	
	//Las CRs de las noches se mantienen salvo que se indique lo contrario para noches de 3+1 (en el caso de buscar hhee, no para cambios)
	if (DesprogramarNoche3mas1 == 1 && BusquedaDeUnCambio==0) {
    	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "n" && CRsNuevas[d][t] >= 4)
    	  	sum(i in Plantilla) (TurneroCalculado[i][d][t]) >= CRsNuevas[d][t]-1;
    	  	
    	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "n" && CRsNuevas[d][t] < 4)
    	  	sum(i in Plantilla) (TurneroCalculado[i][d][t]) >= CRsNuevas[d][t];
  
	} else {
    	forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "n")
    	  	sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]; //No se desprograma ninguna noche	  
	}
           
	//Las CRs para las imag y Oficinas tienen que ser iguales o menores que las del turnero original (así permito desprogramación e impido que se creen nuevas). 
	//Para cuando busco un cambio o hay imag en ese turno (DesprogramarImag=0) esto no lo permito
	if (BusquedaDeUnCambio==0 && DesprogramarImag==1) {
		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "im")
			(sum(i in Plantilla) TurneroCalculado[i][d][t] <= CRsNuevas[d][t]);

		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "it")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) <= CRsNuevas[d][t]);		  

		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "ino")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) <= CRsNuevas[d][t]);
			
		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "o")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) <= CRsNuevas[d][t]);			
	} else {
		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "im")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);

		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "it")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);		  

		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "ino")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);
			
		forall(d in 1..NumDiasTurnero, t in tipoDeTurno : t.Tipo == "o")
			(sum(i in Plantilla) (TurneroCalculado[i][d][t]) == CRsNuevas[d][t]);
	}

	//Dos turnos incompatibles no se pueden asignar (equivale a la restricción de 12h entre dos turnos el mismo día o consecutivos)
      forall(i in Plantilla, tin in turnosIncompatibles, d in 1..NumDiasTurnero-1)
      	(TurneroCalculado[i][d][<tin.t_i>] == 1 => TurneroCalculado[i][d+1][<tin.t_f>] == 0);
     

	//Después de noche o ino tiene que haber 48h excepto si sigue otra noche, ino, v. Tiene en cuenta que Nllm1(im) no se puede hacer. 
	forall(i in Plantilla, d in 1..NumDiasTurnero-3)
		TurneroCalculado[i][d][<"n">] == 1 || TurneroCalculado[i][d][<"ino">] == 1 => 
		((TurneroCalculado[i][d+1][<"l">] == 1) || (TurneroCalculado[i][d+1][<"n">] == 1) || (TurneroCalculado[i][d+1][<"ino">] == 1)) &&
		((TurneroCalculado[i][d+2][<"l">] == 1) || (TurneroCalculado[i][d+2][<"n">] == 1) || (TurneroCalculado[i][d+2][<"ino">] == 1)) &&
		//((TurneroCalculado[i][d+3][<"m1">] == 0) && (TurneroCalculado[i][d+3][<"im">] == 0));
		((sum (t in tipoDeTurno : t.Equiv==1 || t.Equiv == 3) TurneroCalculado[i][d+3][t] ) == 0); //t.Equiv 1, 3 son m1 im
	  

	//Después de doble noche tiene que haber 54h excepto si sigue otra v 
	forall(i in Plantilla, d in 1..NumDiasTurnero-4)
		(TurneroCalculado[i][d][<"n">] == 1 || TurneroCalculado[i][d][<"ino">] == 1) && (TurneroCalculado[i][d+1][<"n">] == 1 || TurneroCalculado[i][d+1][<"ino">] == 1) =>
		(TurneroCalculado[i][d+2][<"l">] == 1) && (TurneroCalculado[i][d+3][<"l">] == 1) &&
		//((TurneroCalculado[i][d+4][<"m1">] == 0) && (TurneroCalculado[i][d+4][<"im">] == 0) && (TurneroCalculado[i][d+4][<"m">] == 0) && (TurneroCalculado[i][d+4][<"o">] == 0));
		((sum (t in tipoDeTurno : t.Equiv>=1 && t.Equiv <= 4) TurneroCalculado[i][d+4][t] ) == 0); // t.Equiv 1 a 4 son m1, im, m, o


		
	//Después de nln tiene que haber 54h excepto si sigue otra v
	forall(i in Plantilla, d in 1..NumDiasTurnero-5, t in tipoDeTurno)
		(TurneroCalculado[i][d][<"n">] == 1 || TurneroCalculado[i][d][<"ino">] == 1) && TurneroCalculado[i][d+1][<"l">] == 1 &&
		(TurneroCalculado[i][d+2][<"n">] == 1 || TurneroCalculado[i][d+2][<"ino">] == 1) =>
		(TurneroCalculado[i][d+3][<"l">] == 1) && (TurneroCalculado[i][d+4][<"l">] == 1) &&
		//((TurneroCalculado[i][d+5][<"m1">] == 0) && (TurneroCalculado[i][d+5][<"im">] == 0) && (TurneroCalculado[i][d+5][<"m">] == 0) && (TurneroCalculado[i][d+5][<"o">] == 0));
		((sum (t in tipoDeTurno : t.Equiv>=1 && t.Equiv <= 4) TurneroCalculado[i][d+5][t] ) == 0); // t.Equiv 1 a 4 son m1, im, m, o


    //Calculo de donde empiezan y acaban los ciclos
    forall(i in Plantilla) Ciclos[i][0]==2 && Ciclos[i][NumDiasTurnero]==2; //día 0 y final siempre 2
    forall(i in Plantilla, d in 1..NumDiasTurnero-1) TurneroCalculado[i][d][<"l">] == 1 => Ciclos[i][d] == 0; //un día libre no puede ser cambio de ciclo
    forall(i in Plantilla, d in 1..NumDiasTurnero-1)
    	(TurneroCalculado[i][d][<"l">] != 1 && TurneroCalculado[i][d+1][<"l">] != 1 => Ciclos[i][d] == 0); //dos turnos consecutivos nunca generan >48h descanso
    
    forall(i in Plantilla, t_r in t_rompeCiclo3D, d in 1..NumDiasTurnero-2)
    	TurneroCalculado[i][d][<t_r.t1>] == 1 && TurneroCalculado[i][d+1][<"l">] == 1 && TurneroCalculado[i][d+2][<t_r.t3>] == 1 =>
    	Ciclos[i][d] == t_r.desc;

    forall(i in Plantilla, t_r in t_rompeCiclo4D, d in 1..NumDiasTurnero-3)
    	TurneroCalculado[i][d][<t_r.t1>] == 1 && TurneroCalculado[i][d+1][<"l">] == 1 && TurneroCalculado[i][d+2][<"l">] == 1 && TurneroCalculado[i][d+3][<t_r.t4>] == 1 =>
    	Ciclos[i][d] == t_r.desc;
   	  

	//Restricciones de ciclo: max numero de dias consecutivos, de mañanas seguidas...  
	forall (i in Plantilla, d in 0..NumDiasTurnero-2, d1 in d+1..NumDiasTurnero-1)
		((Ciclos[i][d]>=1)  && (Ciclos[i][d1] >= 1) && (sum(dx in d+1..d1) (Ciclos[i][dx] >= 1))==1 => 
		(sum(dx2 in d+1..d1, t in tipoDeTurno : t.Tipo != "l") TurneroCalculado[i][dx2][t]) <= 6 && //número de turnos aeronáuticos consecutivos <=6
		(sum(dx2 in d+1..d1, t in tipoDeTurno : t.Tipo == "m" || t.Tipo == "m1" || t.Tipo == "im") TurneroCalculado[i][dx2][t]) <= 5 && //número de mañanas (m, m1, im) seguidas <=5
		(sum(dx2 in d+1..d1, t in tipoDeTurno : t.Tipo == "im" || t.Tipo == "it" || t.Tipo == "ino") TurneroCalculado[i][dx2][t]) <= 2 ); //máximo de dos imaginarias por ciclo

	forall (i in Plantilla, d in 0..NumDiasTurnero-2, d1 in d+1..NumDiasTurnero-1)
		((Ciclos[i][d]>=1 && (Ciclos[i][d1] >= 1) && (sum(dx in d+1..d1) (Ciclos[i][dx] >= 1))==1 && 
		(sum(dx2 in d+1..d1, t in tipoDeTurno : t.Tipo != "l" && t.Tipo != "o") TurneroCalculado[i][dx2][t]) == 6) => Ciclos[i][d1] == 2); //después de 6 días aeronáuticos, descanso >= 54h


//los turnos asignados a un cta son en sentido creciente de tiempo (ordenados temporalmente)

//todos los turnos asignados a un cta tienen que ser diferentes subject to {
  //allDifferent(all(i in R:i%2==1) x[i]);
}

execute DISPLAY {
	for(var i in Plantilla){
		for(var d=1; d<=NumDiasTurnero; d++) {
			for(var t in tipoDeTurno) {
				if (TurneroCalculado[i][d][t] != TurneroActualizado[i][d][t]) {
					write(i," ",t.Tipo," ",d," ",TurneroCalculado[i][d][t]," ",TurneroActualizado[i][d][t],"\n");
  				}
  			}
  		}
  	}				

	//muestra el turnero de la plantilla
	for(var i in Plantilla) {
		write(i," ");
		for(var d=1; d<=NumDiasTurnero; d++) {
			for(var t in tipoDeTurno) {  	
				if (TurneroCalculado[i][d][t] != 0) {
			  		write(t.Tipo, " ");
				}
			}
		}		
		writeln();	
	}	

}

/*main
{
	//cp.param.searchtype = "Restart";
	//cp.param.DynamicProbing="On";


    thisOplModel.generate();

    // Warm start solution
    var sol=new IloOplCPSolution();   
   	for(var i in thisOplModel.Plantilla) {
   	  	for(var d=1; d<=thisOplModel.NumDiasTurnero; d++) {
   			for(var t in thisOplModel.tipoDeTurno){
   				sol.setValue(thisOplModel.TurneroCalculado[i][d][t], thisOplModel.TurneroActualizado[i][d][t]);
   			}
   		}
   	}
    cp.setStartingPoint(sol);

    cp.solve();
    thisOplModel.postProcess();
}

*/


​


    ------------------------------
    [Alex] [Fleischer]
    [EMEA CPLEX Optimization Technical Sales]
    [IBM]
    ------------------------------

    Original Message


  • 7.  RE: User functions in model. Roster model

    Posted Tue July 06, 2021 01:04 PM
    I tried your proposal as a MIP but it didn´t work faster for me. I was trying to implement the warm start for the MIP but I don´t know how to do it. The variable array I want to propose as initial solution is TurneroActualizado

    thisOplModel.generate();
    var vectors = new IloOplCplexVectors();
    vectors.attach(thisOplModel.TurneroCalculado, thisOplModel.TurneroActualizado);
    vectors.setStart(thisOplModel);
    thisOplModel.solve();
    thisOplModel.postProcess();

    could you please help me?

    ------------------------------
    javier rodrigo
    ------------------------------

    Original Message


  • 8.  RE: User functions in model. Roster model

    Posted Wed July 07, 2021 03:04 AM
    Hi,

    for warm start MIP see Warm start through API : addMIPStart within Optimization Simple

        main
    {
    thisOplModel.generate();
    // Warm start the naïve solution
    cplex.addMIPStart(thisOplModel.nbBus,thisOplModel.naiveSolution);
    cplex.solve();
    thisOplModel.postProcess();
    }​


    ------------------------------
    [Alex] [Fleischer]
    [EMEA CPLEX Optimization Technical Sales]
    [IBM]
    ------------------------------

    Original Message


  • 9.  RE: User functions in model. Roster model

    Posted Wed July 07, 2021 07:40 AM
    Thanks. That  worked.

    main
    {
    thisOplModel.generate();
    cplex.addMIPStart(thisOplModel.TurneroCalculado,thisOplModel.TurneroActualizado);
    cplex.solve();
    thisOplModel.postProcess();
    }

    However, it does not make any different to have a warm start or not which is strange because the warm start solution is very similar to the final one. I wonder if it is really taking it into account... I attached a short version of .dat file which finds a solution in 1minute and 25seconds in both cases...

    any thoughts?

    ------------------------------
    javier rodrigo
    ------------------------------

    Original Message


Related Content