Chapter Contents Previous Next
 General Statistics Examples

## Example 8.9: Linear Programming

The two-phase method for linear programming can be used to solve the problem

A routine written in IML to solve this problem follows. The approach appends slack, surplus, and artificial variables to the model where needed. It then solves phase 1 to find a primal feasible solution. If a primal feasible solution exists and is found, the routine then goes on to phase 2 to find an optimal solution, if one exists. The routine is general enough to handle minimizations as well as maximizations.
   /*  Subroutine to solve Linear Programs                      */
/*  names:    names of the decision variables                */
/*  obj:      coefficients of the objective function         */
/*  maxormin: the value 'MAX' or 'MIN', upper or lowercase   */
/*  coef:     coefficients of the constraints                */
/*  rel:      character array of values: '<=' or '>=' or '=' */
/*  rhs:      right-hand side of constraints                 */
/*  activity: returns the optimal value of decision variables*/
/*                                                           */

start linprog( names, obj, maxormin, coef, rel, rhs, activity);

bound=1.0e10;
m=nrow(coef);
n=ncol(coef);

/* Convert to maximization */
if upcase(maxormin)='MIN' then o=-1;
else o=1;

/* Build logical variables */
rev=(rhs<0);
ge=((rel='>=')&^ rev)|((rel='<=')&rev);
eq=(rel='=');
if max(ge)=1 then
do;
sr=I(m);
logicals=-sr[,loc(ge)]||I(m);
artobj=repeat(0,1,ncol(logicals)-m)|(eq+ge);
end;
else do;
logicals=I(m);
artobj=eq;
end;
nl=ncol(logicals);
nv=n+nl+2;


         /* Build coef matrix */
a=((o*obj)||repeat(0,1,nl)||{ -1 0 })//
(repeat(0,1,n)||-artobj||{ 0 -1 })//

/* rhs, lower bounds, and basis */
L=repeat(0,1,nv-2)||-bound||-bound;
basis=nv-(0:nv-1);

/* Phase 1 - primal feasibility */
call lp(rc,x,y,a,b,nv,,l,basis);
print ( { ' ',
'**********Primal infeasible problem************',
' ',
'*********Numerically unstable problem**********',
'*********Singular basis encountered************',
'*******Solution is numerically unstable********',
'***Subroutine could not obtain enough memory***',
'**********Number of iterations exceeded********'
}[rc+1]);
if x[nv] ^=0 then
do;
print '**********Primal infeasible problem************';
stop;
end;
if rc>0 then stop;

/* phase 2 - dual feasibility */
u=repeat(.,1,nv-2)||{ . 0 };
L=repeat(0,1,nv-2)||-bound||0;
call lp(rc,x,y,a,b,nv-1,u,l,basis);

/* Report the solution */
print ( { '*************Solution is optimal***************',
'*********Numerically unstable problem**********',
'**************Unbounded problem****************',
'*******Solution is numerically unstable********',
'*********Singular basis encountered************',
'*******Solution is numerically unstable********',
'***Subroutine could not obtain enough memory***',
'**********Number of iterations exceeded********'
}[rc+1]);
value=o*x  [nv-1];
print ,'Objective Value ' value;
activity= x [1:n] ;
print ,'Decision Variables ' activity[r=names];
lhs=coef*x[1:]:;
dual=y[3:m+2];
print ,'Constraints ' lhs rel rhs dual,
'***********************************************';

finish;

Consider the following product mix example (Hadley 1963). A shop with three machines, A, B, and C, turns out products 1, 2, 3, and 4. Each product must be processed on each of the three machines (for example, lathes, drills, and milling machines). The following table shows the number of hours required by each product on each machine.

 Product Machine 1 2 3 4 A 1.5 1 2.4 1 B 1 5 1 3.5 C 1.5 3 3.5 1

The weekly time available on each of the machines is 2000, 8000, and 5000 hours, respectively. The products contribute 5.24, 7.30, 8.34, and 4.18 to profit, respectively. What mixture of products can be manufactured that maximizes profit? You can solve the problem as follows:

   names={'product 1' 'product 2' 'product 3' 'product 4'};
profit={ 5.24 7.30 8.34 4.18};
tech={ 1.5 1 2.4 1 ,
1 5 1 3.5 ,
1.5 3 3.5 1 };
time={ 2000, 8000, 5000};
rel={ '<=', '<=', '<=' };
run linprog(names,profit,'max',tech,rel,time,products);

The output from this example follows.

 VALUE Objective Value 12737.059

 ACTIVITY Decision Variables product 1 294.11765 product 2 1500 product 3 0 product 4 58.823529

 LHS REL RHS DUAL Constraints 2000 <= 2000 1.9535294 8000 <= 8000 0.2423529 5000 <= 5000 1.3782353

The following example shows how to find the minimum cost flow through a network using linear programming. The arcs are defined by an array of tuples; each tuple names a new arc. The elements in the arc tuples give the names of the tail and head nodes defining the arc. The following data are needed: arcs, cost for a unit of flow across the arcs, nodes, and supply and demand at each node.

A program generates the node-arc incidence matrix and calls the linear program routine for solution:

   arcs={ 'ab' 'bd' 'ad' 'bc' 'ce' 'de' 'ae' };
cost={ 1 2 4 3 3 2 9 };
nodes={ 'a', 'b', 'c', 'd', 'e'};
supdem={ 2, 0, 0, -1, -1 };
rel=repeat('=',nrow(nodes),1);
inode=substr(arcs,1,1);
onode=substr(arcs,2,1);
free n_a_i n_a_o;
do i=1 to ncol(arcs);
n_a_i=n_a_i || (inode[i]=nodes);
n_a_o=n_a_o || (onode[i]=nodes);
end;
n_a=n_a_i - n_a_o;
run linprog(arcs,cost,'min',n_a,rel,supdem,x);

The solution is shown below.

 VALUE Objective Value 8

 ACTIVITY Decision Variables ab 2 bd 2 ad 0 bc 0 ce 0 de 1 ae 0

 LHS REL RHS DUAL Constraints 2 = 2 -2.5 0 = 0 -1.5 0 = 0 -0.5 -1 = -1 -0.5 -1 = -1 -2.5

 Chapter Contents Previous Next Top