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); 
       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 */ 
       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 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:

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 8.9.1.

Output 8.9.1: Product Mix: Optimal Solution

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 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 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.

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 8.9.2.

Output 8.9.2: Minimum Cost Flow: Optimal Solution

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