Originally posted by: Janisch

Hi,

I am completely new to constraint programming. I have already formulated my model into a MILP. Unfortunately, the calculation has taken far too long. Apparently Constraint Programming is much better suited here. Can you help me how I can reformulate my model in CP. Thanks for your help 😊

Model requirements:

It is a sequencing problem (automobile production).

There are 3 variants ("burners", "electric", "gas") of a product. Every variant is produced. For example 50 burners, 40 electric, 50 gas.

There are a fixed number S of stations on a linear assembly line. There are 2 workstations per station. Each product "drives" station by station (station 1 à station 2 à ...) through the assembly.

Each variant claims a processing time at a specific station and a specific workstation. Example:

• burners station 1 Workplace 1 = 118sec

• burners station 1 Workplace 2 = 119sec

• burners station 2 Workplace 1 = 122sec (some processing times exceed the cycle time)

• ……

The cycle time is 120sec.

A station length is also 120sec. However, certain stations may "drift", this means it is may be possible to work longer on a product at the station's two associated workstations. For example, the processing may already take place in the upstream station or the processing may extend into the downstream station!

Only one product may be used at one workstation (not several at the same time!).

The goal is to calculate an optimal mounting sequence that minimizes overload across the entire production line.

I hope the description is sufficient for the beginning :) I really do not know how to model it with CP, with MILP it was quite easy...

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

Hello Janisch,

There are still some ambiguities in the description of the problem.

When you say: "Each variant claims a processing time at a specific station and a specific workstation." what do you mean exactly ? I understand that each product needs to go through all the stations so you have one "activity" for each pair product/station. But what about the 2 workstation of a given station: does the activity uses both, or only one and part of the decisions need to select the workstation that will be used? 

Maybe giving a very small instance with a feasible solution would help.

Philippe

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

By the two work stations I mean, an additional differentiation to each station. Each product enters a station and on the right and left side of the conveyor belt there is a workstation with a worker. This now begins to process the associated process step. The processing times between the workplaces may differ. But both have to be done! I can later upload my MILP wording and an associated dataset. Thanks in advance for the quick reply.

regards

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Now my example:

MILP formulation: I wrote a few comments in the mod- and data file.

Data: Production 3 burners, 3 electric, 3 gas and six stations (export from excel). the second sheet returns the optimal schedule.

x outputs the production sequenc. I also uploaded a screenshot.

Actually I want to count on a total of 210 products and 25 stations

I hope this will help to understand my problem.

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

It helps a little bit but as I do not understand German, the comments are of little help. And I'm also not sure to understand the part around the formulation of drifts and the cost function.

This being said, I think that I could come up with a CP Optimizer model that, hopefully, is quite close to what you want to do. I attach it.

The idea of the model is slightly different from the MIP one in that we do not define intervals for each takt t but for each instance of product (indexed by i in 1..N in my model).

So for each product instance i we have a chain of interval variables station[i][s] representing the "nominal" processing of the product i in station s. These interval variables have a length equal to the station time one (here 120s).

As the stations can work on one product at a time there are sequence variables for each station s (stationSeq[s]) and a noOverlap constraint on these sequences.

The stations work on the different products in the same order (the product cannot bypass each other), so there are some sameSequence constraints between the sequence variables stationSeq. In fact this part of the model is exactly a permutation flowshop problem which is a classical scheduling problem.

Also here as we are in an assembly line context with some takt times, I'm assuming that the "nominal" start time of each work at the station starts at a multiple of the cycle time, which explains the constraint: startOf(station[i][1]) % c == 0.

Now, the products use the different workstations g in 1..G of each station and the duration is not exactly the nominal duration of the stations (120s), and as far as I understand, this is where the drift phenomenon can occur. I create one interval variable work[i][s][g] for the processing of a given product i on workstation g of a station s. Their length depends on the type of product i (known) is given by the data.

If there is no drift time, the intervals work[i][s][g] start at the same date as intervals station[i][s]. The drift time, if any is the distance between the start of station[i][s] and the start of work[i][s][g], this is the meaning of constraint startAtStart(station[i][s], work[i][s][g], drift[i][s][g]). The work stations can work on only one product at a time, so we add noOverlap constraints on the sequence of work of each workstation (workStationSeq). So if a product at a workstation of a station s ends after the nominal end time of stations s for this product, the same work station will have some drift time before starting to work on the next product.

As I do not really understand the objective function, I tried two simple ones: minimize the makespan and minimize the total drift time (given a bound on the makespan, so this seems somehow a bit similar to a limit on the takt times T in the MIP model).

CP Optimizer easily find feasible solutions to the problem, even on problems of the size you mention (210 products and 25 stations) and improves on this initial solution (proofs are provided only for the makespan).

But I'm not sure I modeled the very same problem as yours ...

 

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Hello Phillipe,

Many thanks for your answer. You are very close to what I want to model. However, I am relatively defeated by the model because I do not really understand the CP modeling in OPL yet :)

 

The screenshot drift shows me that a workstation is already starting to work although the upstream workstation has not finished its work yet. I will exclude this.

 

"So for each product instance i we have a chain of interval variables station[i][s] representing the "nominal" processing of the product i in station s. These interval variables have a length equal to the station time one (here 120s)."

Station lengths can also be greater than 120sec if the station allows it. Perhaps it makes sense to introduce a new variable drift range that indicates whether and how many seconds a station can drift.

 

- Is there an indication in which order which Varinate is produced?

- Is it possible to display the variants with a string? Then the evaluation would be clearer.

- Can one sort the ganntchart as shown in the screenshot sequence?

- overload occurs whenever the work can not be completed. This time should be minimized. In practice, for example, a jumper is used to absorb the overload. The model assesses how agile it is positioned for different production programs (for example scenario 1: 100 burners, 20 electrics, 10 gases, scenario 2: 30 burners, 80 electric, 20 gas).

- Would it also be possible to see how much overload occurs at each location or workplace?

 

The model is otherwise very well, what it should. I can not thank you too many times :)

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

You are very close to what I want to model. However, I am relatively defeated by the model because I do not really understand the CP modeling in OPL yet :)

 

Indeed, if you want to use CP Optimizer, you probably should first take a small tour of its modeling (and solving) functionalities. I can suggest this (long) introductory slide deck (https://www.slideshare.net/PhilippeLaborie/introduction-to-cp-optimizer-for-scheduling) and/or this recent survey article (https://ibm.biz/Constraints2018).

 

The screenshot drift shows me that a workstation is already starting to work although the upstream workstation has not finished its work yet. I will exclude this.

 

Indeed, I forgot to add precedence constraints between the consecutive workstations, something like:
forall(i in 1..N, s in 2..S, g in 1..G) { endBeforeStart(work[i][s-1][g],work[i][s][g]);  }   

 

Station lengths can also be greater than 120sec if the station allows it. Perhaps it makes sense to introduce a new variable drift range that indicates whether and how many seconds a station can drift.

 

Possible. In my model the station[] intervals are the nominal intervals at stations, they are synchronized with the take time. And the drift is measured with respect to these ideal start/end times. 

 

- Is there an indication in which order which Varinate is produced?

- Is it possible to display the variants with a string? Then the evaluation would be clearer.

- Can one sort the ganntchart as shown in the screenshot sequence?

 

Just like in CPLEX, you can use the posprocessing in OPL to access the values of the solution and do whatever post-computation you want to display KPIsa about the solution or the solution itself.

 

In CP Optimizer, you can access the start and end values of interval variables (by evaluating expressions startOf(itv), endOf(itv) in the post-processing, by accessing data fields itv.start, itv.end in the OPL script) and sort them or whatever. See for instance the delivered OPL example sched_jobshopflex.

 

- overload occurs whenever the work can not be completed. This time should be minimized. In practice, for example, a jumper is used to absorb the overload. The model assesses how agile it is positioned for different production programs (for example scenario 1: 100 burners, 20 electrics, 10 gases, scenario 2: 30 burners, 80 electric, 20 gas).

- Would it also be possible to see how much overload occurs at each location or workplace?

 

I do not understand the different between overload and the notion of drift time I formulated in the model …

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Thanks for the literature tips. Since I will read in the next few days.

How can I adjust the length of stations in the model? I introduce a new variable "driftlength" [station] [workstation]. I do not want to get this variable in the right place.

 

Now on the subject of overload. A small example:

The assembly line has 4 stations. There is only one variation and we only want to make 1 product. The cycle time is still 120sec.

Processing times [variant] [Station] = [[120], [122], [118], [121]].

Let's say it's about automatic stations. The work can only be done within the station boundaries (length 120sec).

At station 2 and 4 an overload occurs. In the production of only one product a total overload of 3sec.

Now we allow the stations to be able to work across their station boundaries at the work piece (= drift). So there is only 1sec overload. Logically, the end of the assembly line can not be drifted. That's why I want to minimize the overload.

 

I hope this example illustrates the difference between drift and overload.

regards

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

 

How can I adjust the length of stations in the model? I introduce a new variable "driftlength" [station] [workstation]. I do not want to get this variable in the right place.

 

I would create the interval variables stations with a variable length between the nominal length (NominalLength=120) and NominalLength+MaxDrift, something like:

dvar interval station[1..N][s in 1..S] size NominalLength..NominalLength+MaxDrift[s];

 

Then the driftlength is an integer expression:

dexpr int driftlength[i in 1..N][s in 1..S] = lengthOf(station[i][s])-NominalLength;

As a rule of thumb, in CP Optimizer it is better to avoid creating decision variables when you can directly use expressions on existing decision variables (here, you already have the interval variables 'station'). But there are exceptions to this rule.

Now on the subject of overload. A small example:

The assembly line has 4 stations. There is only one variation and we only want to make 1 product. The cycle time is still 120sec.

Processing times [variant] [Station] = [[120], [122], [118], [121]].

Let's say it's about automatic stations. The work can only be done within the station boundaries (length 120sec).

At station 2 and 4 an overload occurs. In the production of only one product a total overload of 3sec.

Now we allow the stations to be able to work across their station boundaries at the work piece (= drift). So there is only 1sec overload. Logically, the end of the assembly line can not be drifted. 

So, if I understand well, the overload is the duration between the end of the station (that can have a drift or not) and the end of the work (if the work ends after the end of the station of course).

So in the model I would formulate it something like maxl(0,endOf(work[i][s][g])-endOf(station[i][s])). 

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Hey Philippe,

nice that I could explain my understanding between drift and overload. I once tried to include the expression in the objective function. I would especially like to define a maximum drift range for each workstation. So far, it is only defined for each station.

Furthermore, I would also like to know how much overload is incurred at each workstation. Can that be modeled?

Unfortunately, the adapted model of you unfortunately does run. But on my PC CP Optimizer could find atfer 10min no solution. Why does it take so long?

I apologize in advance for my inability :)

I attached the adapted model.

regards and a nice weekend

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

Strange. With the model you attached, using the latest version V12.8, a feasible solution is found almost immediately (0.1s) and it is quickly improved until a cost of 55 in less than a couple of seconds. 

Concerning the formulation for the overloads and drift on the work stations, from my understanding, the drift was only defined at the level of the stations as an increase of the nominal station time. So I'm probably still missing something around the definition of drift and overload.

Anyway, in my initial understanding of the problem, so with the curent model, you can just get the overload  at every workstation by considering each term of the cost, so something like:

dexpr int overload[i in 1..N][s in 1..S][g in 1..G] = maxl(0,endOf(work[i][s][g])-endOf(station[i][s]));
minimize sum(i in 1..N, s in 1..S, g in 1..G) overload[i][s][g];

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Hello Philippe,

I have attached my engine log. After 90sec I broke off the calculation. What could it be?

 

- Is it allowed in the model that a workstation can start working earlier if the previous station allows a drift area?

- What may also be missing in the model is the maximum production time of a product. Example: If we have 20 stations, the product is exactly 20 * 120sec in production.

 

regards

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

- The log says that it has found a feasible solution at cost 55. It is just that the engine cannot prove optimality of this solution so it continues until the time limit is reached. But you have a feasible solution.

- Concerning the model of the stations, indeed, the current model is still under the assumption that all stations have a duration equal to 120s and all must start at a date that is a multiple of the takt time. So it should be adapted depending on how exactly the work at the different stations/workstations needs to be synchronized. I suggest you have a look at the presentation/article I mention in order to have a more global view of what you can model with CP Optimizer. It should not be too difficult if you already know CPLEX and MIP because CP Optimizer uses the very same concepts of decision variables, expressions and constraints. It is only that in CP Optimizer you have additional types of decision variables (interval variables, sequence variables), expressions and constraints.

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Is the solution also optimal? There is still a gap of 23%.

Or did I get this wrong.

I'll try to revise the model and then upload it. But if I'm honest, I'm happy that I could set up the model as a MILP. I'm not so good with Cplex yet. I hope you can help me a bit :)

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

CP Optimizer did not prove that the solution is optimal. In fact, in CP Optimizer (and this is maybe a difference compared to MIP), we put a strong focus on (quickly) producing good quality solutions and spend less effort on trying to produce tight lower bounds or optimality proofs.

Philippe

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

In the model file, I have made comments on each line of code. I hope you understand what that means. If something should be unclear, like to ask directly. Maybe then I can answer the question with an example.
Thank you Philipe for your efforts.

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Are my comments in the mod file too short? I am thankful for every help :)

regards

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Hello Philippe,

thanks for your suggestions. I could implement everything so far :)

I understand only the variable maxOverload not yet. Can you give me an example of what I mean by this variable?

I have noticed two errors that I do not know the answer to:

• the following processing times are given (see processing times). If I would now produce a product of each variant and I do not allow drift between jobs (MaxDrift = 0), the calculation should give me an overload of 12sec. But unfortunately that does not come out...

• In the attached image "work times", for some reason, the timing is wrong. everything fits for the first instance (N = 1). But the other instances are no longer in sync. If instance two (N = 2) is produced second, the start times should be 120, 240,360,480 for each of the stations 1-4.

What can not be right there? What do i have to change in the model?

regards

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

> I understand only the variable maxOverload not yet. Can you give me an example of what I mean by this variable?

 

maxOverload is a constant, not a decision variable.

When you write :
dvar interval x size SizeMin..SizeMax;

you define an interval variable whose size is a decision of the problem with a value in the range SizeMin..SizeMax.

The integer expression sizeOf(x) returns the (unknown) size of the interval x that you can use in expressions (like in the cost function) or constraints.

In the model I did:

dvar interval work[i in 1..N][s in 1..S][g in 1..G] size (b[v[i]][s][g]-maxOverload)..b[v[i]][s][g]; 

This means that the size of interval work[i][s][g] is an unknown of the problem between (b[v[i]][s][g]-maxOverload) and b[v[i]][s][g].

And (if I understand well, this is what you want), the overload is therefore b[v[i]][s][g]-sizeOf(work[i][s][g])

 

> In the attached image "work times", for some reason, the timing is wrong. everything fits for the first instance (N = 1). 

> But the other instances are no longer in sync. If instance two (N = 2) is produced second, the start times should be 120, 240,360,480 for each of the stations 1-4.

I think it is because some synchronization of sequences is still missing: synchronize the sequences of the diffrerent workstations of a given station together:

forall(g in 2..G) { 
  sameSequence(stationSeq[1][1],stationSeq[1][g]); 
} 

 

> The following processing times are given (see processing times). If I would now produce a product of each variant and I do not allow drift

> between jobs (MaxDrift = 0), the calculation should give me an overload of 12sec. But unfortunately that does not come out...

I have the feeling that in the current model, the formulation of the first time is kind of redundant, with the "drift" variables an the MaxDrift. In any case, if you force the drift variables to 0, you get an optimal solution with total overload 12:

dvar int drift[i in 1..N][s in 1..S][g in 1..G] in 0..0;//c; 

 

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Hello Philippe,

I hope you read this :)

There are a few things changed about the model requirement:

• It should be possible to define the drift area for each workstation individually, i.e. not for the complete station.
• Example: Only the second station may drift 5 sec. But for some reason, the processing time starts too early at station 4, 5, 6. Station 4 should actually start only at 360sec and station 5 at 480sec and station 6 at 600sec. How can I implement this? See attachment
• At the last station in the assembly line, respective the associated workstations, the processing time must not exceed the assembly line. Example: 3 stations, 2 products, cycle time 120sec. The processing of the first product at the 3rd station (respective the two workstations) can last up to a maximum of 360sec, with the second product up to a maximum of 480sec. --> Perhaps this is feasible with: for all(N,G) endOf (station [i] [S (last station)] [g]) == S + (N-1) * 120

Thank you for your help!

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Update:
I could solve two of the three problems.
However, I still could not solve the problem with the start and end times of the interval variable station:

The drift area (MaxDrift) describes an overlap of two stations. This means that if station s has a maximum drift of 10, then this means that station s is allowed to work into the next station s + 1 for 10s. Conversely, station s + 1 is allowed to start work 10s earlier in station s. Of course, the work should not be superimposed (no parallelism). Depending on the order of production, station s will sometimes work into station s + 1 and sometimes vice versa.

Theoretically, the start times of the stations are fixed if no MaxDrift [S] = 0 is allowed. That is, station 1 starts at 0s, station 2 at 120s, station 3 at 240s, and so on.
However, if Station 2 allows a MaxDrift of 10s (no drift is allowed on all other stations), Station 2 will continue to start at 120s, but may end at 250s, so Station 3 will start at 250s and end at 360s. Suppose Station 2 starts at 120s and ends at 230s, then Station 3 starts at 230s and ends at 360s.


I hope someone has an idea how to model this restriction :)

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: PhilippeLaborie

So it seems that if some pairs of consecutive stations can drift (like stations 2 and 3 in your example), you indeed cannot post the constraints:

startOf(station[i][...]) % 120 == 0;

That forces the station to start "exactely" at times 0,120,240,360,...

You want to say that for some stations i (like the stations 2 and 3 you mention), you can start "around" these multiples of 120, say between n*120-D and n*120+D

So you could do something like:

60-D <= ( 60 + startOf(station[i][...]) ) % 120 <= 60+D

 

 

 

#DecisionOptimization
#OPLusingCPOptimizer

Originally posted by: Janisch

Thank you Phillipe! It works :)

Thanks again for your help :)

 

regards

#DecisionOptimization
#OPLusingCPOptimizer