Choose a Solver

To find the appropriate solver for this problem, consult the Optimization Decision Table. The table asks you to categorize your problem by type of objective function and types of constraints. For this problem, the objective function is linear, and the constraints are linear. The decision table recommends using the linprog solver.

As you see in Problems Handled by Optimization Toolbox Functions or the linprog function reference page, the linprog solver solves problems of the form

The syntax of the linprog solver, as shown in its function reference page, is

[x fval] = linprog(f,A,b,Aeq,beq,lb,ub);

The inputs to the linprog solver are the matrices and vectors in Equation 1.

Combine Variables Into One Vector

There are 16 variables in the equations of Model Description. Put these variables into one vector. The name of the vector of variables is x in Equation 1. Decide on an order, and construct the components of x out of the variables.

The following code constructs the vector using a cell array of names for the variables.

variables = {'I1','I2','HE1','HE2','LE1','LE2','C','BF1',...
    'BF2','HPS','MPS','LPS','P1','P2','PP','EP'};
N = length(variables); 
% create variables for indexing 
for v = 1:N 
   eval([variables{v},' = ', num2str(v),';']); 
end

Executing these commands creates the following named variables in your workspace:

These named variables represent index numbers for the components of x. You do not have to create named variables. The video Optimization Modeling, Part 2: Converting to Solver Form shows how to solve the problem simply using the index numbers of the components of x.

Write Bound Constraints

There are four variables with lower bounds, and six with upper bounds in the equations of Model Description. The lower bounds:

Also, all the variables are positive, which means they have a lower bound of zero.

Create the lower bound vector lb as a vector of 0, then add the four other lower bounds.

lb = zeros(size(variables));
lb([P1,P2,MPS,LPS]) = ...
    [2500,3000,271536,100623];

The variables with upper bounds are:

Create the upper bound vector as a vector of Inf, then add the six upper bounds.

ub = Inf(size(variables));
ub([P1,P2,I1,I2,C,LE2]) = ...
   [6250,9000,192000,244000,62000,142000];

Write Linear Inequality Constraints

There are three linear inequalities in the equations of Model Description:

In order to have the equations in the form A x≤b, put all the variables on the left side of the inequality. All these equations already have that form. Ensure that each inequality is in “less than” form by multiplying through by –1 wherever appropriate:

In your MATLAB^® workspace, create the A matrix as a 3-by-16 zero matrix, corresponding to 3 linear inequalities in 16 variables. Create the b vector with three components.

A = zeros(3,16);
A(1,I1) = 1; A(1,HE1) = -1; b(1) = 132000;
A(2,EP) = -1; A(2,PP) = -1; b(2) = -12000;
A(3,[P1,P2,PP]) = [-1,-1,-1];
b(3) = -24550;

Write Linear Equality Constraints

There are eight linear equations in the equations of Model Description:

In order to have the equations in the form Aeq x=beq, put all the variables on one side of the equation. The equations become:

Now write the Aeq matrix and beq vector corresponding to these equations. In your MATLAB workspace, create the Aeq matrix as an 8-by-16 zero matrix, corresponding to 8 linear equations in 16 variables. Create the beq vector with eight components, all zero.

Aeq = zeros(8,16); beq = zeros(8,1);
Aeq(1,[LE2,HE2,I2]) = [1,1,-1];
Aeq(2,[LE1,LE2,BF2,LPS]) = [1,1,1,-1];
Aeq(3,[I1,I2,BF1,HPS]) = [1,1,1,-1];
Aeq(4,[C,MPS,LPS,HPS]) = [1,1,1,-1];
Aeq(5,[LE1,HE1,C,I1]) = [1,1,1,-1];
Aeq(6,[HE1,HE2,BF1,BF2,MPS]) = [1,1,1,-1,-1];
Aeq(7,[HE1,LE1,C,P1,I1]) = [1267.8,1251.4,192,3413,-1359.8];
Aeq(8,[HE2,LE2,P2,I2]) = [1267.8,1251.4,3413,-1359.8];

Write the Objective

The objective function is

Write this expression as a vector f of multipliers of the x vector:

f = zeros(size(variables));
f([HPS PP EP]) = [0.002614 0.0239 0.009825];

Solve the Problem with linprog

You now have inputs required by the linprog solver. Call the solver and print the outputs in formatted form:

options = optimoptions('linprog','Algorithm','dual-simplex');
[x fval] = linprog(f,A,b,Aeq,beq,lb,ub,options);
for d = 1:N
  fprintf('%12.2f \t%s\n',x(d),variables{d}) 
end
fval

The result:

Optimal solution found.
   136328.74 	I1
   244000.00 	I2
   128159.00 	HE1
   143377.00 	HE2
        0.00 	LE1
   100623.00 	LE2
     8169.74 	C
        0.00 	BF1
        0.00 	BF2
   380328.74 	HPS
   271536.00 	MPS
   100623.00 	LPS
     6250.00 	P1
     7060.71 	P2
    11239.29 	PP
      760.71 	EP

fval =
  1.2703e+03

Examine the Solution

The fval output gives the smallest value of the objective function at any feasible point.

The solution vector x is the point where the objective function has the smallest value. Notice that:

The video Optimization Modeling, Part 2: Converting to Solver Form gives interpretations of these characteristics in terms of the original problem.