# Defining a non-linear optimization-problem

Consider the following non-linear optimisation problem:

WeeklyCapacity = Table[Subscript[maxCapa, t], {t, 6}];
WeeklyDemand = Table[Subscript[weeklyD, i], {i, 10}];
CycleTimes = Table[Subscript[CT, i], {i, 10}];

(*Exemplary data:
WeeklyCapacity = Table[189*7.5, {t, 6}];
WeeklyDemand = Table[RandomInteger[{2000, 5000}], {i, 10}]
CycleTimes = (1/#) & /@ Table[RandomInteger[{40, 125}], {i, 10}]
*)

ProductionProgram =
Table[Subscript[x, i, t], {i, 10}, {t, 6}];

MinVarProductionProgram = Variance/@ ProductionProgram;
ConstraintDemand =
Map[# == 0 &, Total /@ ProductionProgram - WeeklyDemand];
ConstraintCapa = Map[# <= 0 &, Total /@ CapaDemand - WeeklyCapacity];


A company produces i different products and at the end of a week (six working days) has to meet a weekly demand (WeeklyDemand). The production manager knows how many pieces of each product can be produced within one hour (CycleTimes) and is limited to a maximum capacity of hours per day (WeeklyCapacity).

The optimal production program (ProductionProgram) can be calculated as follows:

1. The variance of the production program of each product must be minimised (see MinVarProductionProgram).
2. The production program of each product has to meet its total demand per week (ConstraintDemand)
3. The production program must not exceed the available max. capacity (ConstraintCapa).

MinVarProductionProgram represents a non-linear function. Hence I was thinking about a Lagrangian relaxation approach.

-
Is this homework? – Dr. belisarius May 9 '12 at 23:35
Your objective "function" is a list. Do you want to minimize its sum? The sum of its squares? (That would be nicer because then you'll have a quadratic objective.) For that could use MinStandDeviProductionProgram = Simplify[StandardDeviation /@ ProductionProgram, Assumptions -> Element[Flatten@ProductionProgram, Reals]]; MinStandDeviProductionProgram.MinStandDeviProductionProgram – Daniel Lichtblau May 9 '12 at 23:36
@belisarius: No, it's not homework. I'm trying to level the production in the company I work for. – John May 9 '12 at 23:39
@John Sorry, but the writing sounds like a homework or book exercise. That doesn't mean you will not get answers, but we answer homework-related questions in a more "teach to fish" way. – Dr. belisarius May 9 '12 at 23:39
@belisarius: :D no problem. In fact it's been my intention to write it like that. Sometimes this style is more comprehensive than the abstract way. – John May 9 '12 at 23:44

I'm not sure I fully understand the problem, but maybe this will give you some directions to try.

(*Exemplary data:*)
SeedRandom[11112222333];
WeeklyCapacity = Table[189*7.5, {t, 6}];
WeeklyDemand = Table[RandomInteger[{2000, 5000}], {i, 10}];
CycleTimes = (1/#) & /@ Table[RandomInteger[{40, 125}], {i, 10}];

ProductionProgram = Table[Subscript[x, i, t], {i, 10}, {t, 6}];

ConstraintDemand =
Map[# == 0 &, Total /@ ProductionProgram - WeeklyDemand];
ConstraintCapa = Map[# <= 0 &, Total /@ CapaDemand - WeeklyCapacity];

constraints = Join[ConstraintDemand, ConstraintCapa];
vars = Flatten[ProductionProgram];

MinVarProductionProgram =
Simplify[Variance /@ ProductionProgram,
Assumptions -> Element[vars, Reals]];


Here I am not sure whether it is the total or total-of-squares (or something else entirely) to be minimized.

Timing[{min, vals} =
FindMinimum[{Total[MinVarProductionProgram], constraints}, vars]]

Out[294]= {0.1, {-8.53712*10^-11, {Subscript[x, 1, 1] -> 389.5,
Subscript[x, 1, 2] -> 389.5, Subscript[x, 1, 3] -> 389.5,
Subscript[x, 1, 4] -> 389.5, Subscript[x, 1, 5] -> 389.5,
Subscript[x, 1, 6] -> 389.5, Subscript[x, 2, 1] -> 549.167,
Subscript[x, 2, 2] -> 549.167, Subscript[x, 2, 3] -> 549.167,
Subscript[x, 2, 4] -> 549.167, Subscript[x, 2, 5] -> 549.167,
Subscript[x, 2, 6] -> 549.167, Subscript[x, 3, 1] -> 703.667,
Subscript[x, 3, 2] -> 703.667, Subscript[x, 3, 3] -> 703.667,
Subscript[x, 3, 4] -> 703.667, Subscript[x, 3, 5] -> 703.667,
Subscript[x, 3, 6] -> 703.667, Subscript[x, 4, 1] -> 495.167,
Subscript[x, 4, 2] -> 495.167, Subscript[x, 4, 3] -> 495.167,
Subscript[x, 4, 4] -> 495.167, Subscript[x, 4, 5] -> 495.167,
Subscript[x, 4, 6] -> 495.167, Subscript[x, 5, 1] -> 759.833,
Subscript[x, 5, 2] -> 759.833, Subscript[x, 5, 3] -> 759.833,
Subscript[x, 5, 4] -> 759.833, Subscript[x, 5, 5] -> 759.833,
Subscript[x, 5, 6] -> 759.833, Subscript[x, 6, 1] -> 764.167,
Subscript[x, 6, 2] -> 764.167, Subscript[x, 6, 3] -> 764.167,
Subscript[x, 6, 4] -> 764.167, Subscript[x, 6, 5] -> 764.167,
Subscript[x, 6, 6] -> 764.167, Subscript[x, 7, 1] -> 637.333,
Subscript[x, 7, 2] -> 637.333, Subscript[x, 7, 3] -> 637.333,
Subscript[x, 7, 4] -> 637.333, Subscript[x, 7, 5] -> 637.333,
Subscript[x, 7, 6] -> 637.333, Subscript[x, 8, 1] -> 476.5,
Subscript[x, 8, 2] -> 476.5, Subscript[x, 8, 3] -> 476.5,
Subscript[x, 8, 4] -> 476.5, Subscript[x, 8, 5] -> 476.5,
Subscript[x, 8, 6] -> 476.5, Subscript[x, 9, 1] -> 666.833,
Subscript[x, 9, 2] -> 666.833, Subscript[x, 9, 3] -> 666.833,
Subscript[x, 9, 4] -> 666.833, Subscript[x, 9, 5] -> 666.833,
Subscript[x, 9, 6] -> 666.833, Subscript[x, 10, 1] -> 511.,
Subscript[x, 10, 2] -> 511., Subscript[x, 10, 3] -> 511.,
Subscript[x, 10, 4] -> 511., Subscript[x, 10, 5] -> 511.,
Subscript[x, 10, 6] -> 511.}}}


For sum of squares of variances:

Timing[{min2, vals2} =
FindMinimum[{MinVarProductionProgram.MinVarProductionProgram,
constraints}, vars]]

Out[296]= {2.28, {2.51657*10^8, {Subscript[x, 1, 1] -> 478.496,
Subscript[x, 1, 2] -> 373.593, Subscript[x, 1, 3] -> 373.593,
Subscript[x, 1, 4] -> 373.593, Subscript[x, 1, 5] -> 373.593,
Subscript[x, 1, 6] -> 364.134, Subscript[x, 2, 1] -> 686.262,
Subscript[x, 2, 2] -> 523.509, Subscript[x, 2, 3] -> 523.509,
Subscript[x, 2, 4] -> 523.509, Subscript[x, 2, 5] -> 523.509,
Subscript[x, 2, 6] -> 514.701, Subscript[x, 3, 1] -> 877.119,
Subscript[x, 3, 2] -> 668.976, Subscript[x, 3, 3] -> 668.976,
Subscript[x, 3, 4] -> 668.976, Subscript[x, 3, 5] -> 668.976,
Subscript[x, 3, 6] -> 668.979, Subscript[x, 4, 1] -> 346.338,
Subscript[x, 4, 2] -> 552.386, Subscript[x, 4, 3] -> 552.386,
Subscript[x, 4, 4] -> 552.386, Subscript[x, 4, 5] -> 552.386,
Subscript[x, 4, 6] -> 415.116, Subscript[x, 5, 1] -> 861.299,
Subscript[x, 5, 2] -> 743.336, Subscript[x, 5, 3] -> 743.336,
Subscript[x, 5, 4] -> 743.336, Subscript[x, 5, 5] -> 743.336,
Subscript[x, 5, 6] -> 724.357, Subscript[x, 6, 1] -> 937.424,
Subscript[x, 6, 2] -> 727.414, Subscript[x, 6, 3] -> 727.414,
Subscript[x, 6, 4] -> 727.414, Subscript[x, 6, 5] -> 727.414,
Subscript[x, 6, 6] -> 737.922, Subscript[x, 7, 1] -> 775.12,
Subscript[x, 7, 2] -> 611.552, Subscript[x, 7, 3] -> 611.552,
Subscript[x, 7, 4] -> 611.552, Subscript[x, 7, 5] -> 611.552,
Subscript[x, 7, 6] -> 602.671, Subscript[x, 8, 1] -> 612.536,
Subscript[x, 8, 2] -> 451.053, Subscript[x, 8, 3] -> 451.053,
Subscript[x, 8, 4] -> 451.053, Subscript[x, 8, 5] -> 451.053,
Subscript[x, 8, 6] -> 442.253, Subscript[x, 9, 1] -> 766.477,
Subscript[x, 9, 2] -> 650.705, Subscript[x, 9, 3] -> 650.705,
Subscript[x, 9, 4] -> 650.705, Subscript[x, 9, 5] -> 650.705,
Subscript[x, 9, 6] -> 631.703, Subscript[x, 10, 1] -> 523.304,
Subscript[x, 10, 2] -> 523.034, Subscript[x, 10, 3] -> 523.034,
Subscript[x, 10, 4] -> 523.034, Subscript[x, 10, 5] -> 523.034,
Subscript[x, 10, 6] -> 450.56}}}


Hope this gives some ideas for how to proceed.

--- edit ---

Since the objective is nonlinear Mathematica only has NMinimize to try to enforce integrality of variables. Here is the altered code for this situation. I start by rounding the result from FindMinimum, to be used as initial variable ranges for NMinimize.

In[35]:= Timing[{min2, vals2} =
FindMinimum[{MinVarProductionProgram.MinVarProductionProgram,
constraints}, vars];]

Out[35]= {2.16, Null}

In[39]:= firstGuess = Round[vars /. vals2];
delta = 50;
ranges = Transpose[{vars, firstGuess - delta, firstGuess + delta}];


I use these ranges in NMinimize.

Timing[{min3, vals3} =
NMinimize[{MinVarProductionProgram.MinVarProductionProgram,
Append[constraints, Element[vars, Integers]]}, ranges,
MaxIterations -> 1000]]

During evaluation of In[42]:= NMinimize::cvmit: Failed to converge to the requested accuracy or precision within 1000 iterations. >>

Out[42]= {153.86, {5.87444, {Subscript[x, 1, 1] -> 389,
Subscript[x, 1, 2] -> 390, Subscript[x, 1, 3] -> 390,
Subscript[x, 1, 4] -> 389, Subscript[x, 1, 5] -> 389,
Subscript[x, 1, 6] -> 390, Subscript[x, 2, 1] -> 548,
Subscript[x, 2, 2] -> 550, Subscript[x, 2, 3] -> 547,
Subscript[x, 2, 4] -> 550, Subscript[x, 2, 5] -> 550,
Subscript[x, 2, 6] -> 550, Subscript[x, 3, 1] -> 704,
Subscript[x, 3, 2] -> 704, Subscript[x, 3, 3] -> 704,
Subscript[x, 3, 4] -> 705, Subscript[x, 3, 5] -> 703,
Subscript[x, 3, 6] -> 702, Subscript[x, 4, 1] -> 495,
Subscript[x, 4, 2] -> 495, Subscript[x, 4, 3] -> 495,
Subscript[x, 4, 4] -> 496, Subscript[x, 4, 5] -> 495,
Subscript[x, 4, 6] -> 495, Subscript[x, 5, 1] -> 759,
Subscript[x, 5, 2] -> 759, Subscript[x, 5, 3] -> 761,
Subscript[x, 5, 4] -> 760, Subscript[x, 5, 5] -> 760,
Subscript[x, 5, 6] -> 760, Subscript[x, 6, 1] -> 764,
Subscript[x, 6, 2] -> 763, Subscript[x, 6, 3] -> 764,
Subscript[x, 6, 4] -> 765, Subscript[x, 6, 5] -> 765,
Subscript[x, 6, 6] -> 764, Subscript[x, 7, 1] -> 638,
Subscript[x, 7, 2] -> 638, Subscript[x, 7, 3] -> 636,
Subscript[x, 7, 4] -> 638, Subscript[x, 7, 5] -> 637,
Subscript[x, 7, 6] -> 637, Subscript[x, 8, 1] -> 477,
Subscript[x, 8, 2] -> 476, Subscript[x, 8, 3] -> 477,
Subscript[x, 8, 4] -> 476, Subscript[x, 8, 5] -> 476,
Subscript[x, 8, 6] -> 477, Subscript[x, 9, 1] -> 666,
Subscript[x, 9, 2] -> 666, Subscript[x, 9, 3] -> 667,
Subscript[x, 9, 4] -> 667, Subscript[x, 9, 5] -> 668,
Subscript[x, 9, 6] -> 667, Subscript[x, 10, 1] -> 511,
Subscript[x, 10, 2] -> 511, Subscript[x, 10, 3] -> 511,
Subscript[x, 10, 4] -> 511, Subscript[x, 10, 5] -> 511,
Subscript[x, 10, 6] -> 511}}}


As the message indicates, possibly one could do better. Notice though that the min is now considerably lower than what we had from FindMinimum, so progress has been made in the globval optimization effort. And of course we can keep going. This time I'll narrow the start range lengths.

nextGuess = vars /. vals3;
delta2 = 10;
ranges2 = Transpose[{vars, nextGuess - delta2, nextGuess + delta2}];

Timing[{min4, vals4} =
NMinimize[{MinVarProductionProgram.MinVarProductionProgram,
Append[constraints, Element[vars, Integers]]}, ranges2,
MaxIterations -> 1000]]

Out[66]= {135.86, {0.461111, {Subscript[x, 1, 1] -> 389,
Subscript[x, 1, 2] -> 389, Subscript[x, 1, 3] -> 390,
Subscript[x, 1, 4] -> 390, Subscript[x, 1, 5] -> 390,
Subscript[x, 1, 6] -> 389, Subscript[x, 2, 1] -> 549,
Subscript[x, 2, 2] -> 549, Subscript[x, 2, 3] -> 550,
Subscript[x, 2, 4] -> 549, Subscript[x, 2, 5] -> 549,
Subscript[x, 2, 6] -> 549, Subscript[x, 3, 1] -> 704,
Subscript[x, 3, 2] -> 704, Subscript[x, 3, 3] -> 704,
Subscript[x, 3, 4] -> 703, Subscript[x, 3, 5] -> 703,
Subscript[x, 3, 6] -> 704, Subscript[x, 4, 1] -> 496,
Subscript[x, 4, 2] -> 495, Subscript[x, 4, 3] -> 495,
Subscript[x, 4, 4] -> 495, Subscript[x, 4, 5] -> 495,
Subscript[x, 4, 6] -> 495, Subscript[x, 5, 1] -> 760,
Subscript[x, 5, 2] -> 760, Subscript[x, 5, 3] -> 760,
Subscript[x, 5, 4] -> 759, Subscript[x, 5, 5] -> 760,
Subscript[x, 5, 6] -> 760, Subscript[x, 6, 1] -> 764,
Subscript[x, 6, 2] -> 765, Subscript[x, 6, 3] -> 764,
Subscript[x, 6, 4] -> 764, Subscript[x, 6, 5] -> 764,
Subscript[x, 6, 6] -> 764, Subscript[x, 7, 1] -> 637,
Subscript[x, 7, 2] -> 638, Subscript[x, 7, 3] -> 638,
Subscript[x, 7, 4] -> 637, Subscript[x, 7, 5] -> 637,
Subscript[x, 7, 6] -> 637, Subscript[x, 8, 1] -> 477,
Subscript[x, 8, 2] -> 476, Subscript[x, 8, 3] -> 477,
Subscript[x, 8, 4] -> 476, Subscript[x, 8, 5] -> 477,
Subscript[x, 8, 6] -> 476, Subscript[x, 9, 1] -> 666,
Subscript[x, 9, 2] -> 667, Subscript[x, 9, 3] -> 667,
Subscript[x, 9, 4] -> 667, Subscript[x, 9, 5] -> 667,
Subscript[x, 9, 6] -> 667, Subscript[x, 10, 1] -> 511,
Subscript[x, 10, 2] -> 511, Subscript[x, 10, 3] -> 511,
Subscript[x, 10, 4] -> 511, Subscript[x, 10, 5] -> 511,
Subscript[x, 10, 6] -> 511}}}


Seems to be stabilizing.

--- end edit ---

-
This looks pretty good! The only thing which is missing, is that all variables must to be element of Integer. – John May 10 '12 at 15:54
Lucky thing, because it seems the above doesn't find good optima anyway. Will provide an edit. – Daniel Lichtblau May 10 '12 at 16:14