GA or Surrogateopt doesn't find solutions | global nonlinear integer optimization in factory planning

Hi,

I have some problems with my optimization model.

I'm trying to minimize the costs resulting from the needed capacity to produce products in a factory. So for example I'm deciding on how much machines, workers, work shift, overtimes, on which systems the products should be produced and so on. I need to find the global minimum and have nonlinear constraints as well as integer objective variables.

I'm using the global optimization toolbox and tried the genetic algorithm(it always just finds a local minimum which doesn't fulfill the constraints) and the surrogateopt(which doesn't find a solution but my initial point solution).

With surrogateopt I can see, that his wild guess objectivevariables are not fulfilling my constraints. How can I modify the solver maybe to a branch-and-bound or branch-and-cut learning?

So I think either I have to many objective Variables with to many dependencies, I'm calling the objective function in a wrong way or my constraints are to constricting/wrong.

Could you maybe try to help me with my problem?

I will try to show you my code in a very summarized form so that it does not get too complex with its 4500 rows.

Main script:

objVar = 18;

% objVar for each period(t) includes Decisionvariables [(1)BeMi Sys 1, (2)BeMi Sys2, (3)Ausstattung1, (4)Ausstattung2, (5)LagerA, (6)LagerB,

% (7)FertigwarenlagerA, (8)FertigwarenlagerB, (9)Überstunden Sys1, (10)Schichtmodell, (11)Überstunden Sys2,

% (12)Fremdverbage A, (13)Fremdverbage B, (14)BeMisMinSys1 (15)BeMisMinSys1 (16)neueBeMisSys1 (17)neueBeMisSys2

% (18)BouleanChangeShift]

numVar = objVar * 14; %at the moment T=14, but needs to be 26 up to 52

Intcon = [1:numVar]; %all my objVariables are Integer since I e.g. can just buy a whole Machine and not a half

lbvar18Periode1 = [1 1 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0];

lbvar18 = [1 1 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0];

lb = [lbvar18Periode1, repmat(lbvar18,1,VPerioden-1)];

ubvar18Periode1 = [15 10 3 3 6 6 10 10 20 3 20 0 0 0 0 7 3 0];

ubvar18 = [15 10 3 3 6 6 10 10 20 3 20 0 0 7 3 7 3 1];

ub = [ubvar18Periode1, repmat(ubvar18,1,VPerioden-1)];

options = optimoptions('surrogateopt',...

'Display', 'iter',...

'PlotFcn', {'optimplotfvalconstr', 'optimplotfval', 'optimplotx', 'surrogateoptplot'},...

'CheckpointFile', 'checkfile.mat',...

'MaxFunctionEvaluations', 5000,... % Default is max(200,50*nvar), For debugging I scaled down the Evaluations, A bad solution is fine to evaluate the code at the moment

'UseParallel', false,... % Parallel Computing doesn't work at the moment. I donn't know why...

'InitialPoints', Startpoint... % I'm giving some initialpoints, so that there are solutions (obviously not optimal) the solver can learn from.

);

objFunction = @(OptimVar)objconstr(OptimVar,Nachfrage, solution_p,l1,l2,t1_p,t2_p,FL1,FL2);

[solution,objectiveValue, exitflag] = surrogateopt(objFunction, lb, ub, Intcon, options);

I only use nonlinear constraints.

Objconstr function for surrogateopt:

function f = objconstr(OptimVar,Nachfrage, solution_p,l1,l2,t1_p,t2_p,FL1,FL2)

global M MAKosten BeMiKosten Ausstattungskosten Ueberstunde Materiallagerkosten Fertiglagerkosten...

K_Fremdvergabe Kosten_Material tstop objVar VPerioden Schichtsys ZVLAnlage MEZAnlage MEZSchicht grosseZahlI;

%%First some constraints like:

% ConstraintType1: Worktime at 5 workdays can not overcome 120h

j = 10;

t=0;

while t <= VPerioden - 1

%Amount of WorkShifts(1to3)*8h/shift*5days + OvertimeInWeek <= 120

f.Ineq(i) = OptimVar(j+t*objVar)*8*5 + OptimVar((j-1)+t*objVar) - 120;

i = i + 1;

t = t + 1;

end

%ConstraintType2: For constraint minimum time to hold same shift-model

t=0;

while t < VPerioden-1

% |AmountOfShifts(t=2)-AmountOfShifts(t=2)| <= BouleanChangeShift(t=2)*M=10;

% So my BouleanChangeShift shows if there was a change in shift-model (changing shift Model is expensive)

f.Ineq(i) = abs(OptimVar(10+(t+1)*objVar) - OptimVar(10+t*objVar)) - (OptimVar(18+(t+1)*objVar)*M);

i = i + 1;

t = t + 1;

end

t=1;

while t < VPerioden-1

if t <= (VPerioden-MEZSchicht) %(VPerioden=14; MEZSchicht=8)

% Explanation: When the Shift Model changes: the next 8 periods the shift model has to be the same size to fulfill be se sum as 8 and so to fulfill the constraint

% I think this constraint is very hard to fulfill with random values for all these objvalues ?! What do you think?

% BouleanChangeShift(t=2)* ( AmountOfShifts(t=2)/AmountOfShifts(t=3) + AmountOfShifts(t=3)/AmountOfShifts(t=4) + AmountOfShifts(t=4)/AmountOfShifts(t=5)

% + AmountOfShifts(t=5)/AmountOfShifts(t=6) + AmountOfShifts(t=6)/AmountOfShifts(t=7) + AmountOfShifts(t=8)/AmountOfShifts(t=2)

% <= BouleanChangeShift(t=2)*MEZSchicht

f.Ineq(i) = OptimVar(18+t*objVar) * (OptimVar(10+t*objVar)/OptimVar(10+(t+1)*objVar)+OptimVar(10+(t+1)*objVar)/OptimVar(10+(t+2)*objVar)+...

OptimVar(10+(t+2)*objVar)/OptimVar(10+(t+3)*objVar)+OptimVar(10+(t+3)*objVar)/OptimVar(10+(t+4)*objVar)+...

OptimVar(10+(t+4)*objVar)/OptimVar(10+(t+5)*objVar)+OptimVar(10+(t+5)*objVar)/OptimVar(10+(t+6)*objVar)+...

OptimVar(10+(t+6)*objVar)/OptimVar(10+(t+7)*objVar)+OptimVar(10+(t+7)*objVar)/OptimVar(10+t*objVar)) - ...

(OptimVar(18+t*objVar) * MEZSchicht);

i = i + 1; %% I need equality of the shift model so I need the workarround with A <= B and A >= B:

f.Ineq(i) = -(OptimVar(18+t*objVar) * (OptimVar(10+t*objVar)/OptimVar(10+(t+1)*objVar)+OptimVar(10+(t+1)*objVar)/OptimVar(10+(t+2)*objVar)+...

OptimVar(10+(t+2)*objVar)/OptimVar(10+(t+3)*objVar)+OptimVar(10+(t+3)*objVar)/OptimVar(10+(t+4)*objVar)+...

OptimVar(10+(t+4)*objVar)/OptimVar(10+(t+5)*objVar)+OptimVar(10+(t+5)*objVar)/OptimVar(10+(t+6)*objVar)+...

OptimVar(10+(t+6)*objVar)/OptimVar(10+(t+7)*objVar)+OptimVar(10+(t+7)*objVar)/OptimVar(10+t*objVar)) - ...

(OptimVar(18+t*objVar) * MEZSchicht));

i = i + 1;

elseif t == (VPerioden-(MEZSchicht-1))

(...)

elseif t == (VPerioden-(MEZSchicht-6))

f.Ineq(i) = OptimVar(18+t*objVar) * (OptimVar(10+t*objVar)/OptimVar(10+(t+1)*objVar)+OptimVar(10+(t+1)*objVar)/OptimVar(10+t*objVar)) - ...

(OptimVar(18+t*objVar) * MEZSchicht);

i = i + 1;

f.Ineq(i) = -(OptimVar(18+t*objVar) * (OptimVar(10+t*objVar)/OptimVar(10+(t+1)*objVar)+OptimVar(10+(t+1)*objVar)/OptimVar(10+t*objVar)) - ...

(OptimVar(18+t*objVar) * MEZSchicht));

i = i + 1;

end

t = t + 1;

end

%% Now my objfunction to minimize he overall cost of all periods like:

% Kostenberechnung für jede Periode im Betrachtungszeitraum

Anlagenausstattung1 = zeros(2, VPerioden);

t=1;

while t <= VPerioden

switch OptimVar(3+(t-1)*objVar) %see how System1 is equiped and take into account the corresponding cost rate of vector "Ausstattungskosten"

case 1 %Ausstattungsvariante 1 --> Kostensatz Ausstattungskosten2

Anlagenausstattung1(:,t) = [OptimVar(3+(t-1)*objVar), Ausstattungskosten(2)];

case 2 %Ausstattungsvariante 2 --> Kostensatz Ausstattungskosten3

Anlagenausstattung1(:,t) = [OptimVar(3+(t-1)*objVar), Ausstattungskosten(3)];

case 3 %Ausstattungsvariante 3 --> Kostensatz Ausstattungskosten4

Anlagenausstattung1(:,t) = [OptimVar(3+(t-1)*objVar), Ausstattungskosten(4)];

end

t=t+1;

end

%initialize the costvectors

cost_Mitarbeiter = zeros(1,VPerioden);

cost_UeberstundenMA = zeros(1,VPerioden);

cost_Ressourcen = zeros(1,VPerioden);

cost_Periode = zeros(1,VPerioden);

%calculate the cost per period

while t <= VPerioden

% Mitarbeiterkosten = AnzAnlagen*2*Schichtsystem*Betr.KOSTENSATZ_MA

cost_Mitarbeiter(t) = Mitarbeiter_Anzahl(t)* MAKosten(2)*40;

% ÜberstundenkostenMA = AnzÜberstundenAnlageA*AnzAnlagen*2MA/Anlage*KOSTENSATZ_Überstunde

cost_UeberstundenMA(t) = OptimVar(9+(t-1)*objVar)* OptimVar(1+(t-1)*objVar)* 2* Ueberstunde(1)+ OptimVar(11+(t-1)*objVar)* OptimVar(2+(t-1)*objVar)* 2* Ueberstunde(1);

% Betriebskosten der Kapazitäten-Betriebsmittel, Materiallagereinheiten, Fertigwarenlagereinheiten

cost_Ressourcen(t) = (OptimVar(1+(t-1)*objVar)* Anlagenausstattung1(2,t)* BeMiKosten(2)* 8 *OptimVar(10+(t-1)*objVar)* 5)...

+(OptimVar(2+(t-1)*objVar)* Anlagenausstattung2(2,t)* BeMiKosten(2)* 8 *OptimVar(10+(t-1)*objVar)* 5) ...

+(OptimVar(5+(t-1)*objVar)+ OptimVar(6+(t-1)*objVar)) * Materiallagerkosten(2)*50 ...

+(OptimVar(7+(t-1)*objVar)+ OptimVar(8+(t-1)*objVar)) * Fertiglagerkosten(2)*50; %+ OptimVar(12)* Fremdvergabe %Kosten Schichtmodell durch Anzahl notwendige Mitarbeiter

cost_Periode(t) = cost_Mitarbeiter(t)+ cost_UeberstundenMA(t)+ cost_Ressourcen(t);

t=t+1;

end

%Call functions to plan the production of products.

[Produktion1, Produktion2, Nachfrage, l1, l2] = Produktionsplanung(Nachfrage, OptimVar,l1,l2,t1_p,t2_p,FL1,FL2);

[~, ~, Produktion1, Produktion2] = zeit_Verschiebung(Produktion1, Produktion2);

[MehrBelastung,x2,x7,Produktion12,Produktion22, Nachfrage,l1,l2] = ort_Verschiebung(OptimVar,Produktion1,Produktion2,Nachfrage,l1,l2,FL1,FL2);

[Auftrag_fremd, Nachfrage] = Fremdvergabe(Nachfrage,l1,l2,FL1,FL2, Produktion12, Produktion22, OptimVar);

Produktion12 = [Produktion12; Auftrag_fremd(Auftrag_fremd(:,1) == 1,1:7)];

Produktion22 = [Produktion22; Auftrag_fremd(Auftrag_fremd(:,1) == 2,1:7)];

[y3,y4] = Verspaetung(Produktion12, Produktion22);

y1 = y3 + y4;

[s2,~] = size(Auftrag_fremd);

Kosten_Fremd = s2 * K_Fremdvergabe;

Nachfrage = Nachfrage(Nachfrage(:,1)>0,1:4);

f1 = find(Nachfrage(:,3) < tstop);

c1 = f1./f1;

s1 = sum(c1);

k1 = s1 * highcostrate;

%And this is the final cost value and the objectivevalue for the optimization

f.Fval = sum(cost_Periode)+ sum(inv_Periode)+ sum(cost_Material)+ y1 +k1+ Kosten_Fremd;

At the Moment I'm considering 14 Periods. The Problem has a size of 203 constraints an 252 objectiveVariables. As I mentioned I have to consider 26 to 52 Periods, so the problime will have to be scaled up.

By calling the surrogateopt I'm displaying the actual vector of the objVariables in every iteration and I can see, that they don't fullfil my constraints (except my initialpoints) even not after 10,000 iterations.

I also tried a lot of options for the GA (with the corresponding modified code) from default to personal, but I never got a good result.

Is it permissible to call the objectivefunction like this(function handle)? When I run my code row by row the GA doesn't call the objective function for every population.

Has anybody a good hint, how I can solve my problem? E.G. use different solver, change script,...?

How can I use non-objective variables in the constraints, that need to change depending on decision variables? see Example "BouleanChangeShift" This Boulean depends in the decision of the other objective-Variables "AmountOfShifts(t=2)-AmountOfShifts(t=2)". It woulb be great if this could be a calculated variable and then is out of the decision space of the solver. So the problem would be smaller:

% |AmountOfShifts(t=2)-AmountOfShifts(t=2)| <= BouleanChangeShift(t=2)*M=10;

% So my BouleanChangeShift shows if there was a change in shift-model (changing shift Model is expensive)

f.Ineq(i) = abs(OptimVar(10+(t+1)*objVar) - OptimVar(10+t*objVar)) - (OptimVar(18+(t+1)*objVar)*M);

I would be extremely grateful if you could help me, since I'm stuck in this problem for a month now.

Thanks a lot in advance!

Answer 1

[Alan Weiss]

You asked a lot of questions, but it seems to me that what you really want to know is whether there are any feasible points for your problem. That is a hard question to answer in general. Here are some steps you can take, as described in Converged to an Infeasible Point:

• Remove your objective function and nonlinear constraints and see whether there is any point satisfying the integer and linear constraints by calling intlinprog. If there are no feasible points, then you know that your formulation is too restrictive.
• If there are feasible points to the linear + integer constrained problem, remove the integer constraints and replace the nonlinear constraints. Call fmincon to see whether there are any solutions to the problem without integer constraints. This can take some experimentation because you might need to try several initial points.
• If you get feasible points with both of these relaxations, then you might obtain some values of points to try as initial points for ga or surrogateopt. But there are no guarantees.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Niklas Rochow]

Hi Alan,

Thank you for your advice and sorry for that many questions. I think its because I'm stuck in the problem for so long time now.

First according to my questions, I will try to deepen the question:

My objectivefunction has a lot of calculations in it using the objectivevariables the solver is trying and is also calling functions to at the end get the objectivevalue e.g. under the comment

% Betriebskosten der Kapazitäten-Betriebsmittel, Materiallagereinheiten, Fertigwarenlagereinheiten

and

%Call functions to plan the production of products.

At the End of my objectivefunction the objectuvevalue is calculated by

f.Fval = sum(cost_Periode)+ sum(inv_Periode)+ sum(cost_Material)+ y1 +k1+ Kosten_Fremd;

So in difference to the examples of the MathWorks Documentation there is no direct link of the objectivevalue and the objectivevariables in the objectivefunction. Mathworks examples always show cases like

objval = x1*cost1 + x2*cost2

My guess is, that this results in a bad solution of the solver. Am I right and do I have to change my objectivefunction? But I don't know how to change it.

Now according to your solution/advice:

I tried sizing down my problem by leaving just 14 objectivevariables all of the same kind within 1 difference of upper (ub=7) and lower bound(lb=6)

1. removing all my constraints -->Result: surrogateopt finds solution (my initialpoint is best)
2. removing nonlinear constraints -->Result: surrogateopt finds solution (my initialpoint is best)
3. trying with all constraints -->Result: surrogateopt finds solution (my initialpoint is best)

When I'm sizing up my problem by giving greater decision space (ub=14; lb=2) he doesn't find solution but my initialpoint.

This leads me to my guess.

Thank you for your help!