Example 9.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*/
   /*                                                           */
   *ods trace output;
   proc iml;
   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);
      adj=(-1*rev)+^ rev;
      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 */
      *ods trace output;
      proc iml;
      a=((o*obj)||repeat(0,1,nl)||{ -1 0 })//
        (repeat(0,1,n)||-artobj||{ 0 -1 })//
        ((adj#coef)||logicals||repeat(0,m,2));

         /* rhs, lower bounds, and basis */
      b={0,0}//(adj#rhs);
      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:n];
      dual=y[3:m+2];
      print ,'Constraints ' lhs rel rhs dual,
             '***********************************************';

   finish;

Consider the following product mix example (Hadley; 1962). 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:

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 results from this example are shown in Output 9.9.1.

The following example shows how to find the minimum cost flow through a network by 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 that define 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.

The following 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 in Output 9.9.2.