Set Up a Linear Program, Problem-Based

Convert a Problem to Solver Form

This example shows how to convert a linear problem from mathematical form into Optimization Toolbox™ solver syntax using the problem-based approach.

The variables and expressions in the problem represent a model of operating a chemical plant, from an example in Edgar and Himmelblau [1]. There are two videos that describe the problem.

The remainder of this example is concerned solely with transforming the problem to solver syntax. The example closely follows the video Optimization Modeling, Part 2: Problem-Based Solution of a Mathematical Model.

Model Description

The video Mathematical Modeling with Optimization, Part 1 suggests that one way to convert a problem into mathematical form is to:

  1. Get an overall idea of the problem

  2. Identify the goal (maximizing or minimizing something)

  3. Identify (name) variables

  4. Identify constraints

  5. Determine which variables you can control

  6. Specify all quantities in mathematical notation

  7. Check the model for completeness and correctness

For the meaning of the variables in this section, see the video Mathematical Modeling with Optimization, Part 1.

The optimization problem is to minimize the objective function, subject to all the other expressions as constraints.

The objective function is:

0.002614 HPS + 0.0239 PP + 0.009825 EP.

The constraints are:

2500P16250
I1192,000
C62,000
I1 - HE1132,000
I1 = LE1 + HE1 + C
1359.8 I1 = 1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P1
3000P29000
I2244,000
LE2142,000
I2 = LE2 + HE2
1359.8 I2 = 1267.8 HE2 + 1251.4 LE2 + 3413 P2
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
LPS = LE1 + LE2 + BF2
MPS = HE1 + HE2 + BF1 - BF2
P1 + P2 + PP24,550
EP + PP12,000
MPS271,536
LPS100,623
All variables are positive.

First Solution Method: Create Optimization Variable for Each Problem Variable

The first solution method involves creating an optimization variable for each problem variable. As you create the variables, include their bounds.

P1 = optimvar('P1','LowerBound',2500,'UpperBound',6250);
P2 = optimvar('P2','LowerBound',3000,'UpperBound',9000);
I1 = optimvar('I1','LowerBound',0,'UpperBound',192000);
I2 = optimvar('I2','LowerBound',0,'UpperBound',244000);
C = optimvar('C','LowerBound',0,'UpperBound',62000);
LE1 = optimvar('LE1','LowerBound',0);
LE2 = optimvar('LE2','LowerBound',0,'UpperBound',142000);
HE1 = optimvar('HE1','LowerBound',0);
HE2 = optimvar('HE2','LowerBound',0);
HPS = optimvar('HPS','LowerBound',0);
MPS = optimvar('MPS','LowerBound',271536);
LPS = optimvar('LPS','LowerBound',100623);
BF1 = optimvar('BF1','LowerBound',0);
BF2 = optimvar('BF2','LowerBound',0);
EP = optimvar('EP','LowerBound',0);
PP = optimvar('PP','LowerBound',0);

Create Problem and Objective

Create an optimization problem container. Include the objective function in the problem.

linprob = optimproblem('Objective',0.002614*HPS + 0.0239*PP + 0.009825*EP);

Create and Include Linear Constraints

There are three linear inequalities in the problem expressions:

I1 - HE1132,000
EP + PP12,000
P1 + P2 + PP24,550.

Create these inequality constraints and include them in the problem.

linprob.Constraints.cons1 = I1 - HE1 <= 132000;
linprob.Constraints.cons2 = EP + PP >= 12000;
linprob.Constraints.cons3 = P1 + P2 + PP >= 24550;

There are eight linear equalities:

I2 = LE2 + HE2
LPS = LE1 + LE2 + BF2
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
I1 = LE1 + HE1 + C
MPS = HE1 + HE2 + BF1 - BF2
1359.8 I1 = 1267.8 HE1 + 1251.4 LE1 + 192 C + 3413 P1
1359.8 I2 = 1267.8 HE2 + 1251.4 LE2 + 3413 P2.

Include these constraints as well.

linprob.Constraints.econs1 = LE2 + HE2 == I2;
linprob.Constraints.econs2 = LE1 + LE2 + BF2 == LPS;
linprob.Constraints.econs3 = I1 + I2 + BF1 == HPS;
linprob.Constraints.econs4 = C + MPS + LPS == HPS;
linprob.Constraints.econs5 = LE1 + HE1 + C == I1;
linprob.Constraints.econs6 = HE1 + HE2 + BF1 == BF2 + MPS;
linprob.Constraints.econs7 = 1267.8*HE1 + 1251.4*LE1 + 192*C + 3413*P1 == 1359.8*I1;
linprob.Constraints.econs8 = 1267.8*HE2 + 1251.4*LE2 + 3413*P2 == 1359.8*I2;

Solve Problem

The problem formulation is complete. Solve the problem using solve.

linsol = solve(linprob);
Optimal solution found.

Examine Solution

Evaluate the objective function. (You could have asked for this value when you called solve.)

evaluate(linprob.Objective,linsol)
ans =

   1.2703e+03

The lowest-cost method of operating the plant costs $1,207.30.

Examine the solution variable values.

tbl = struct2table(linsol)
tbl =

  1×16 table

    BF1    BF2      C         EP         HE1           HE2           HPS            I1           I2       LE1       LE2           LPS           MPS         P1       P2       PP  
    ___    ___    ______    ______    __________    __________    __________    __________    ________    ___    __________    __________    __________    ____    ______    _____

    0      0      8169.7    760.71    1.2816e+05    1.4338e+05    3.8033e+05    1.3633e+05    2.44e+05    0      1.0062e+05    1.0062e+05    2.7154e+05    6250    7060.7    11239

This table is too wide to see easily. Stack the variables to get them to a vertical orientation.

vars = {'P1','P2','I1','I2','C','LE1','LE2','HE1','HE2',...
'HPS','MPS','LPS','BF1','BF2','EP','PP'};
outputvars = stack(tbl,vars,'NewDataVariableName','Amt','IndexVariableName','Var')
outputvars =

  16×2 table

    Var       Amt    
    ___    __________

    P1          6250
    P2        7060.7
    I1    1.3633e+05
    I2      2.44e+05
    C        8169.7
    LE1             0
    LE2    1.0062e+05
    HE1    1.2816e+05
    HE2    1.4338e+05
    HPS    3.8033e+05
    MPS    2.7154e+05
    LPS    1.0062e+05
    BF1             0
    BF2             0
    EP        760.71
    PP         11239
  • BF1, BF2, and LE1 are 0, their lower bounds.

  • I2 is 244,000, its upper bound.

  • The nonzero components of the objective function (cost) are

    • HPS380,328.74

    • PP11,239.29

    • EP760.71

The video Optimization Modeling, Part 2: Problem-Based Solution of a Mathematical Model gives interpretations of these characteristics in terms of the original problem.

Second Solution Method: Create One Optimization Variable and Indices

Alternatively, you can solve the problem using just one optimization variable that has indices with the names of the problem variables. This method enables you to give a lower bound of zero to all problem variables at once.

vars = {'P1','P2','I1','I2','C','LE1','LE2','HE1','HE2',...
    'HPS','MPS','LPS','BF1','BF2','EP','PP'};
x = optimvar('x',vars,'LowerBound',0);

Set Variable Bounds

Include the bounds on the variables using dot notation.

x('P1').LowerBound = 2500;
x('P2').LowerBound = 3000;
x('MPS').LowerBound = 271536;
x('LPS').LowerBound = 100623;
x('P1').UpperBound = 6250;
x('P2').UpperBound = 9000;
x('I1').UpperBound = 192000;
x('I2').UpperBound = 244000;
x('C').UpperBound = 62000;
x('LE2').UpperBound = 142000;

Create Problem, Linear Constraints, and Solution

The remainder of the problem setup is similar to the setup using separate variables. The difference is that, instead of addressing a variable by its name, such as P1, you address it using its index, x('P1').

Create the problem object, include the linear constraints, and solve the problem.

linprob = optimproblem('Objective',0.002614*x('HPS') + 0.0239*x('PP') + 0.009825*x('EP'));

linprob.Constraints.cons1 = x('I1') - x('HE1') <= 132000;
linprob.Constraints.cons2 = x('EP') + x('PP') >= 12000;
linprob.Constraints.cons3 = x('P1') + x('P2') + x('PP') >= 24550;

linprob.Constraints.econs1 = x('LE2') + x('HE2') == x('I2');
linprob.Constraints.econs2 = x('LE1') + x('LE2') + x('BF2') == x('LPS');
linprob.Constraints.econs3 = x('I1') + x('I2') + x('BF1') == x('HPS');
linprob.Constraints.econs4 = x('C') + x('MPS') + x('LPS') == x('HPS');
linprob.Constraints.econs5 = x('LE1') + x('HE1') + x('C') == x('I1');
linprob.Constraints.econs6 = x('HE1') + x('HE2') + x('BF1') == x('BF2') + x('MPS');
linprob.Constraints.econs7 = 1267.8*x('HE1') + 1251.4*x('LE1') + 192*x('C') + 3413*x('P1') == 1359.8*x('I1');
linprob.Constraints.econs8 = 1267.8*x('HE2') + 1251.4*x('LE2') + 3413*x('P2') == 1359.8*x('I2');

[linsol,fval] = solve(linprob);
Optimal solution found.

Examine Indexed Solution

Examine the solution as a vertical table.

tbl = table(vars',linsol.x')
tbl =

  16×2 table

    Var1        Var2   
    _____    __________

    'P1'           6250
    'P2'         7060.7
    'I1'     1.3633e+05
    'I2'       2.44e+05
    'C'          8169.7
    'LE1'             0
    'LE2'    1.0062e+05
    'HE1'    1.2816e+05
    'HE2'    1.4338e+05
    'HPS'    3.8033e+05
    'MPS'    2.7154e+05
    'LPS'    1.0062e+05
    'BF1'             0
    'BF2'             0
    'EP'         760.71
    'PP'          11239

Bibliography

[1] Edgar, Thomas F., and David M. Himmelblau. Optimization of Chemical Processes. McGraw-Hill, New York, 1988.

Related Topics