Originally posted by: qtbgo

Hi, All.

I have the follewing CP code:

forall( i in vessels, t in T) t >= s[i] && t < e[i] => r3[i][t] >= r_min[i] && r3[i][t] <= r_max[i];   forall( i in vessels, t in T) t < s[i] || t >= e[i] => r3[i][t] == 0;

Where s, e, r3 are decision variables. r_min, r_max are data.
What I want is when t is in s_i to e_i , then r3_it should be in a range, otherwise r3_it should be 0.

The above code works, but very inefficient when the problem is large. I test and find that the second constraint takes the longest time.

Does anyone have more efficient solution to this situation?
Thanks in advance.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: GGR

Hi

the solver report should give you the number of constraint you add the model. In your sample it comes card(vessels)*card(T). I figure out vessels is an enumeration of vessels. But what is the meaning of T.

Just a remark, if T is a time (which seems quite natural) and of big dimensions, you could consider an interval based scheduling model.

Hope that helps
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

GGR, vessels is 1..40.
T is, for example,1..500. It's time.

I am trying an interval based scheduling model, but not finished yet.

Now, this non interval based scheduling model works, but inefficient for large problems.

I am considering ways to improve its performance.

First, I saw somewhere OPL sentence like
search { try endtry}, but cannot find any documents on it in online help.

second, I think there may be other ways to model this situation which is more efficient, so I request help here.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: GGR

Hi

the search { try endtry} is no more available in OPL. All that stuff is now replaced by the notion of search phase.

Basically, in your example, you possibly may to choose deciding first on vessel allocation or time allocation: for this you will use a search phase.

About the scheduling model, the main point of interest in to avoid a complete enumeration of the time. The search of CP Optimizer works like that (unless very specific cases as proof of optimality). So you must avoid any expression and constraint set that enumerate the time (in practice you must not have to write forall(t in T) It is at this condition you can hope having better result using an interval variable based model.

Cheers
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Thanks GGR.

I have figured out an interval based model, but it is even more inefficient. So, I still work on the first model. The three main constraints are:

forall( i in vessels, t in T) t >= s[i] && t < e[i] => r3[i][t] >= r_min[i] && r3[i][t] <= r_max[i]; forall( i in vessels, t in T) t < s[i] || t >= e[i] => r3[i][t] == 0;   forall(i,j in vessels : i<j) (s[j] < e[i] && s[j] >= s[i]) || (s[i] < e[j] && s[i] >= s[j]) => !((b[j] < b [i] + l[i] && b[j] >= b[i]) ||(b[i] < b [j] + l[j] && b[i] >= b[j]) ); 
//noOverlap simultaneously, b is also decision variable.

According to your experience, which one is not good? And any suggestion to improve it?
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: GGR

Hi

I think you must adapt your whole problem to a complete interval based model. that is,

Your two first constraints clearly define a cumulative constraint:

int R_Len_max = ...;   

int N = ...; 

int R_min[1..N] = ...; 

int R_max[1..N] = ...; 

int R_Len_vessels[1..N] = ...;   dvar interval R_vessels[i in 1..N] size R_Len_vessels[i];   cumulFunction Len_Vessels = sum (i in 1..N) pulse(R_vessels[i], R_min[i], R_max[i]);   dexpr 

int R_expr[i in 1..N] = heigthAtStart(Len_Vessels, R_Vessels[i])   subject to 
{ R_Len_vessels <= R_Len_max; 
}

Then your disjunction is quite complicated; you can have it simpler by writing:

forall(i,j in vessels : i<j) (e[j] <= s[i]) || (e[i] >= s[j]) || (b[j] + l[j] <= b[i]) || (b[i] + l[i] <= b[j])

Last point, Your models looks like 2 Dimensional placement on a the temporal and a B cumulative dimension. that is the [b, b + l)'s are interval variables. If I am right, and if there is a cumulative limitation on the B dimension, it is a good idea to have a cumulative expression on it.

Hope that helps
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Yes , it is in nature a 2 dimensional placement on a temporal and a space dimension.
My english is poor, so I cannot express it clearly.

My scheduling based model is as follows(Main parts)

using CP; 

int n = ... ; 
//number of vessels 

int L = ... ; 
// number of 10-meter berth segments 

int H = ... ; 
// number of 1-hour time periods 

int Q = ... ;  
// number of available cranes   range T = 0..H-1; range vessels = 1..n;   

int l[vessels] = ... ; 
// length of vessel 

int ETA[vessels] = ... ; 
// expected time of arrival 

int EST[vessels] = ... ; 
// earliest starting time 

int EFT[vessels] = ... ; 
// expected finishing time 

int r_min[vessels] = ... ; 
// minimum number of cranes to assign 

int r_max[vessels] = ... ; 
// maximum number of cranes to assign   dvar interval vTask[i in vessels][t in T] optional in EST[i]..H size 1; 
// time dimension. A vessel's handling time is divided into a serial of 1 hour segments. Each segment may require different number of Quays.  dvar interval vt[i in vessels] in 0..H ;  
//for span vTask    dvar interval vb[i in vessels] in 0..L  size l[i];  
//space dimension. The berthing place of vessels   dexpr 

int  s[i in vessels] = min(t in T) startOf(vTask[i][t], H); 
//start of handling of a vessel dexpr 

int  e[i in vessels] = max(t in T) endOf(vTask[i][t],0); 
//end of handling of a vessel       cumulFunction totalQuaysUsage = sum(i in vessels, t in T) pulse(vTask[i][t],r_min[i], r_max[i]);   
//dexpr int individualQuayUsage[i in vessels][t in T] = heightAtStart(vTask[i][t], totalQuaysUsage);//   subject to
{   totalQuaysUsage <= Q ; forall(i in vessels, t in T) presenceOf(vTask[i][t]) => startOf(vTask[i][t]) == t; 
//fix vTask[i][t] 
//NoOverlap on time space diagram forall(i in vessels ) span(vt[i], all(t in T)vTask[i][t]); forall( i,j in vessels: i!=j ) overlapLength(vb[i], vb[j]) !=0 => overlapLength(vt[i], vt[j]) ==0 ; 
}

This code can only solve 3 vessels problem. It cannot give even a feasible solution when there are 10 vessels which can be solved in MIP formulation by Cplex.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

In fact, the problem I am addressing is described in the attached paper.

I want to test if CP method can surpass the heuristics in the paper.

Note of the moderator : paper attachment has been removed.

Edited by: ilog_opt_admin on Jan 30, 2012 3:40 AM
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: SystemAdmin

Hello,
I understand from the paper that the number of quay cranes assigned to a vessel can change every time units. In your model you create an interval variable per time-unit in 0..horizon per vessel. It would be better to split the interval variables representing the vessels into sub-intervals of size 1 given a duration min and max you can compute on the vessel. This will result in less interval variables and they can be more tightly constrained with the vessel time intervals. In particular the number of interval variables will not depend on the schedule horizon.

Here is the corresponding model where I also introduced some minimal handling of the total number of quay cranes assigned to a vessel.

The unit intervals for each vessel are constrained as a chain of optional intervals with endAtStart constraints. If the minimal and maximal duration of vessels are dMin and dMax (computed from the vessel crane capacity demand and the max and min number of quay cranes r_max, r_min), the first dMin unit intervals are constrained to be present and the other ones are such that the ith one implies the i-1 one.

I also introduced a redundant cumulative constraint for the berth length and a search phase to first fix the time-related variables (vut) before the exact positioning on the berth (vb).

I think this model should work better.

using CP;   

int n = ...; 
// number of vessels 

int L = ...; 
// number of 10-meter berth segments 

int H = ...; 
// number of 1-hour time periods 

int Q = ...; 
// number of available cranes   range vessels = 1..n; 

int l    [vessels] = ... ; 
// length of vessel 

int ETA  [vessels] = ... ; 
// expected time of arrival 

int EST  [vessels] = ... ; 
// earliest starting time 

int EFT  [vessels] = ... ; 
// expected finishing time 

int r_min[vessels] = ... ; 
// minimum number of cranes to assign 

int r_max[vessels] = ... ; 
// maximum number of cranes to assign 

int m    [vessels] = ... ; 
// crane capacity demand of vessel given as a number of QC-hours   range T = 0..H-1;   

int dMin[v in vessels] = (m[v]+r_max[v]-1) div r_max[v]; 
// (m div r_max) rounded up 

int dMax[v in vessels] = (m[v]+r_min[v]-1) div r_min[v]; 
// (m div r_min) rounded up   tuple VesselUnit 
{ 

int vessel; 

int number; 
} 
{VesselUnit
} vUnits = 
{ <v, n> | v in vessels, n in 1..dMax[v] 
};   
// Time dimension dvar interval vut[u in vUnits] optional in EST[u.vessel]..H size 1; 
// A vessel's handling time is divided into a serial of 1 hour segments. Each segment may require different number of Quays.  dvar interval vt[v in vessels] in EST[v]..H size dMin[v]..dMax[v]; 
// Space dimension dvar interval vb[v in vessels] in 0..L size l[v]; 
// Space dimension. The berthing place of vessels 
// Quay crane usage 
// Alternative model that factorizes r_min on the vessel interval cumulFunction totalQuaysUsage = sum(v in vessels) pulse(vt[v], r_min[v]) + sum(u in vUnits) pulse(vut[u], 0, r_max[u.vessel]-r_min[u.vessel]); dexpr 

int individualQuayUsage[u in vUnits] = r_min[u.vessel]+heightAtStart(vut[u], totalQuaysUsage);   cumulFunction berthLengthUsage = sum(v in vessels) pulse(vt[v], l[v]); execute 
{ var f = cp.factory; cp.setSearchPhases(f.searchPhase(vut), f.searchPhase(vb)); 
}   minimize max(v in vessels) endOf(vt[v]); subject to 
{ 
// Quay crane number forall(v in vessels) 
{ span(vt[v], all(u in vUnits: u.vessel==v) vut[u]); m[v] <= sum(u in vUnits: u.vessel==v) individualQuayUsage[u] <= m[v]+r_min[v]; forall(u in vUnits: u.vessel==v) 
{ 

if (1<u.number) 
{ endAtStart(vut[<v,u.number-1>], vut[u]); 
} 

if (u.number <= dMin[v]) 
{ presenceOf(vut[u]); 
} 

else 
{ presenceOf(vut[u])=>presenceOf(vut[<v,u.number-1>]); 
} 
} 
} totalQuaysUsage <= Q ; 
// NoOverlap on time space diagram forall( i,j in vessels: i!=j ) 
{ overlapLength(vt[i], vt[j]) * overlapLength(vb[i], vb[j]) == 0; 
} 
// Redundant cumulative constraint berthLengthUsage <= L; 
}

Philippe
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Think you so much, Philippe.

I combined your code into my final model to test. It's quick now, but cannot give correct solution.

For now, It seems that

dexpr 

int individualQuayUsage[u in vUnits] = r_min[u.vessel]+heightAtStart(vut[u], totalQuaysUsage);

should be

dexpr 

int individualQuayUsage[u in vUnits] = presenceOf(vut[u])? r_min[u.vessel]: 0 +heightAtStart(vut[u], totalQuaysUsage);

But this will lead to no solution in my model.

In the attachment, there is a MIP and a CP model for the same problem. You can test them(data is also supplied).
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: SystemAdmin

Yes, you are right the formulation of the individualQuayUsage expression is not correct in my model: if the interval vut[u] is absent, the quay crane usage at the time unit should be 0 and not r_min[u.vessel].
Note that your formulation is also wrong because of operators priorities:

presenceOf(vut\[u\])? r_min\[u.vessel\]: 0 +heightAtStart(vut\[u\], totalQuaysUsage);

is interprested as:

presenceOf(vut\[u\])? r_min\[u.vessel\]: ( 0 +heightAtStart(vut\[u\], totalQuaysUsage) );

whereas you need:

(presenceOf(vut\[u\])? r_min\[u.vessel\]: 0) + heightAtStart(vut\[u\], totalQuaysUsage);

or

(presenceOf(vut[u])*r_min[u.vessel]) + heightAtStart(vut[u], totalQuaysUsage);

Note that you can also avoid that by using the simpler model:

cumulFunction totalQuaysUsage = sum(u in vUnits) pulse(vut[u], r_min[u.vessel], r_max[u.vessel]); dexpr 

int individualQuayUsage[u in vUnits] = heightAtStart(vut[u], totalQuaysUsage);

But this model looses the factorization of the r_min quantity over the whole duration of the vessel. I don't know which model is better. The best of to try both.

On both fixed models, CP Optimizer now finds a solution with same cost as the CPLEX one (6.5).

Note that CP Optimizer will have difficulties to prove optimality even for very small instances for this type of problem. Here, you should mainly use it to compute good solutions.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Thank you, Philippe.

It works now. Now I report results to you.
I computed a big instance in the attachment.
a TS heuristics can solve it within 4 minutes with an objective of 254.
Whereas our CP method cannot even find a feasible solution within at least 15 minutes.

I think maybe this problem is not suitable to solve using CP, Or is there any way to improve it further?
By the way, It seems that

int dMax[v in vessels] = (m[v]+r_min[v]-1) div r_min[v] ; 
// (m div r_min) rounded up

should be

int dMax[v in vessels] = ftoi(minl(ceil(m[v]* (1 + beta*L) /(r_min[v]^alpha)), H)); 
// considering berth deviation factor and interference exponent when there are multiple quays

But it still cannot solve big problem.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: SystemAdmin

It is probably because the horizon is a bit tight for the computation of an initial feasible solution.
You can try increasing the H value in the data. Unlike in the MIP model, this will have no effect on the size of the CP Optimizer model. Then, the objective function (lateness and tardiness costs) should in any case minimize the end time of the schedule.
I tried doing that in the model for the large data you attached (using a very big H=1680 instead of H=168), then CP Optimizer finds an initial solution (of quite poor quality) much more easily and within the 4mn time-limit it reach a solution with cost 240. And that solution indeed satisfies a schedule horizon of 168. I'm using V12.4 and the new version for dMax that you proposed. I attach the model but it should be the same as the one you are running.
If the tight horizon H is a strong constraint in the real problem and not an artefact to limit the size of the MIP model, you can try solving the problem in two steps: first step to find a solution with minimal makespan (and hopefully makespan <= H) with a version of the model that minimizes only the makespan (starting from a larger horizon) followed by a second steps that minimizes the real objective. The solution of the first step can be re-injected as starting point to the second step (see notion of Starting Point in CP Optimizer).
Note that the computation of dMax has an impact on the model size so any upper bound in the definition of dMax will help to reduce the size of the model.
Another potential source of difficulty to produce an initial solution with the present model could be the search phase that first fixes the time variables before the spatial berth positions. This phase is very useful to produce compact solutions on the time axis using a cumulative relaxation (berthLengthUsage) of the limited berth length but it may lead to difficulties in finding feasible berth positions. One thing you can do if you encounter this type of problem is also to exchange the order of the search phase (fix intervals vb before the vut) in an initial search step to produce a feasible solution.
But, at least on the large instance you attached, it is not clear whether those initial steps are necessary.
Philippe
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Good results.
But I cannot reproduce your result. I changed H to 1680, but still cannot find a feasible solution.

Maybe this is because I'm using version 12.2. I'm now applying for an acdemic version 12.4 to see if it works.

By the way, could you delete the paper I attached before(becauseof the copyright)?

Thanks very much.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Dear Philipple,
I have got Cplex 12.3 academic version (cannot get v12.4). I test the big instance with your CP.mod . I set H to 1068. But it still cannot find a feasible solution. I don't know why.
Do you change something that I don't know?

Any hints?
thanks.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: SystemAdmin

I attach the version of the model and data I'm using, just in case they are different from the ones you are using. I checked and it easily finds a solution both with V12.3 and V12.4.
If the projects turn out to be the same, could you send the exact build version of your V12.3 (you can find it in the windows that pops u when you click on: Help > About IBM ILOG CPLEX Optimization Studio) and the begining of the search log?
Here is the begining of the search log I get using V12.3:

!

---

! Minimization problem - 1171 variables, 6007 constraints, 2 phases
! Preprocessing : 807 extractables eliminated
! TimeMode = ElapsedTime
! TimeLimit = 240
! Initial process time : 0,04s (0,00s extraction + 0,04s propagation)
! . Log search space : 11874,5 (before), 11874,2 (after)
! . Memory usage : 6,9 MB (before), 11,5 MB (after)
! Using parallel search with 2 workers.
!

---

! Best Branches Non-fixed W Branch decision
1000 584 1 on vut#423
1000 583 2 on vut#996
2000 18 1 on vb(13)
2000 44 2 on vut#660
* 7544,4 2017 1,34s 1 -
7544,4 3000 23 2 on vut#426
7544,4 3000 593 1 on vut#670
7544,4 4000 23 2 on vut#715
7544,4 4000 61 1 on vut#962
* 7467,8 4075 1,81s 1 -
7467,8 5000 23 2 on vut#955
7467,8 5000 635 1 on vut#1002
Philippe
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

Sorry, I made a mistake. Now it works as yours.

Thank you very much.
#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: qtbgo

just to check my question has been answered
#DecisionOptimization
#OPLusingCPOptimizer