Details Examples (omodd01)

/***************************************************************/
/*                                                             */
/*          S A S   S A M P L E   L I B R A R Y                */
/*                                                             */
/*    NAME: omodd01                                            */
/*   TITLE: Details Examples (omodd01)                         */
/* PRODUCT: OR                                                 */
/*  SYSTEM: ALL                                                */
/*    KEYS: OR                                                 */
/*   PROCS: OPTMODEL                                           */
/*    DATA:                                                    */
/*                                                             */
/* SUPPORT:                             UPDATE:                */
/*     REF:                                                    */
/*    MISC: Examples from the Details section of the OPTMODEL  */
/*          chapter of Mathematical Programming.               */
/*                                                             */
/***************************************************************/

/*  Named Parameters  */
data coeff;
   input c_xx c_x c_y c_xy c_yy;
   datalines;
1 -1  -2 -1 1
;
proc optmodel;
   var x, y;
   number c_xx, c_x, c_y, c_xy, c_yy;
   read data coeff into c_xx c_x c_y c_xy c_yy;
   min z=c_xx*x**2 + c_x*x + c_y*y + c_xy*x*y + c_yy*y**2;
   solve;
quit;

/*  Indexing  */
proc optmodel;
   number p{1..3};
   p[1]=5;
   p[2]=7;
   p[3]=9;
   put p[*]=;
quit;

/*  Indexing in an INIT Clause  */
proc optmodel;
   number p{i in 1..3} init 3 + 2*i;
   put (sum{i in 1..3} p[i]);
quit;

/*  Indexing with a SUM Operator  */
proc optmodel;
   number n init 100000;
   var x{1..n};
   min z = sum{i in 1..n}(x[i] - log(i))**2;
   solve;
quit;

/*  Parameters  */
proc optmodel;
   number pi=4*atan(1);
   number r;
   number circum=2*pi*r;
   r=1;
   put circum;       /* prints 6.2831853072 */
   r=2;
   put circum;       /* prints 12.566370614 */
quit;

/*  Expressions  */
proc optmodel;
   var x init 0.5 >= 0 <= 1;
   var y init (0.5 >= 0) <= 1;
   print x y;
   expand;
quit;

proc optmodel;
   number original{i in 1..8} = sin(i);
   number sorted{i in 1..8} init original[i];
   call sortn(of sorted[*]);
   print original sorted;
quit;

/*  Index Sets  */
proc optmodel;
   put (sum{i in 1..10} i**2);
quit;

/*  Nil Values  */
proc optmodel;
   num a{1..2};
   a[1] = 1;
   put (a[1] = _nil_)=; /* false, defined */
   put (a[2] = _nil_)=; /* true, undefined */
   put ((if a[1] then _nil_ else 2) = a[1])= /* false */;
   put ((if a[1] then _nil_ else 2) = a[2])= /* true */;
   a[1] = a[2]; /* now a[1] is undefined */
   put (a[1] = _nil_)=; /* true, undefined */

   /* array initialization list is evaluated at most once */
   num b{1..3} init [1 2 3 [*] 7];
   put b[*];     /* evaluates the initialization list */
   b[1] = _nil_; /* reset b[1] */
   put b[*];     /* b[1] replaced by the default value */

   var x init (rand('UNIFORM')) >= 0 <= 10;
   put x=;       /* initialized to a random value */
   x = _nil_;    /* discard the current value */
   put x=;       /* initialized to a new random value */
   put x.ub=;    /* the declared bound */
   x.ub = 1;     /* override the declared bound */
   put x.ub=;    /* the overridden bound */
   x.ub = _nil_; /* reverse the override */
   put x.ub=;    /* the declared bound again */

   /* nil values print as empty cells */
   print {i in 1..3, j in 1..3}(if i+j-1=3 then _nil_ else 1/(i+j-1));
quit;

/*  AND Aggregation Expression  */
proc optmodel;
   put (and{i in 1..5} i < 10); /* returns 1 */
   put (and{i in 1..5} i NE 3); /* returns 0 */
quit;

/*  CARD Function  */
proc optmodel;
   put (card(1..3));
quit;

/*  CROSS Expression  */
proc optmodel;
   set s1 = 1..2;
   set<string> s2 = {'a', 'b'};
   set<number, string> s3=s1 cross s2;
   put 's3 is ' s3;
   set<number, string, number> s4 = s3 cross 4..5;
   put 's4 is ' s4;
quit;

/*  DIFF Expression  */
proc optmodel;
   put ({1,3} diff {2,3}); /* outputs {1} */
quit;

/*  Equality Expression  */
proc optmodel;
   put (/a b/ = /a b/);  /* true */
   put (/a b/ = /a/);    /* false */
   put ({} = /a b/);     /* false */
   put (1..3 = {1,3,2}); /* true, ignores order */
quit;

/*  IN Expression  */
proc optmodel;
   set s = 1..10;
   put (5 in s);         /* outputs 1 */
   put (-1 not in s);    /* outputs 1 */
   set<num, str> t = {<1,'a'>, <2,'b'>, <2,'c'>};
   put (<2, 'b'> in t);   /* outputs 1 */
   put (<1, 'b'> in t);  /* outputs 0 */
quit;

/*  Index Set Expression  */
proc optmodel;
   put ({i in 1..5 : i NE 3}); /* outputs {1,2,4,5} */
quit;

/*  INTER Expression  */
proc optmodel;
   put ({1,3} inter {2,3}); /* outputs {3} */
quit;

/*  INTER Aggregation Expression  */
proc optmodel;
   put (inter{i in 1..3} i..i+3); /* outputs {3,4} */
quit;

/*  MAX Aggregation Expression  */
proc optmodel;
   put (max{i in 2..5} 1/i);
   put (max{i in 2..5} repeat('*',i-1));
quit;

/*  MIN aggregation Expression  */
proc optmodel;
   put (min{i in 2..5} 1/i);
   put (min{i in 2..5} repeat('*',i-1));
quit;

/*  OR Aggregation Expression  */
proc optmodel;
   put (or{i in 1..5} i = 2); /* returns 1 */
   put (or{i in 1..5} i = 7); /* returns 0 */
quit;

/*  PROD Aggregation Expression  */
proc optmodel;
   number n = 5;
   put (prod{i in 1..n} i); /* outputs 120 */
quit;

/*  RANGE Expression  */
proc optmodel;
   put (10..30 by 7); /* outputs {10,17,24} */
quit;

/*  Set Constructor Expression  */
proc optmodel;
   put ({1,2,3,2}); /* outputs {1,2,3} */
quit;

/*  Display of Tuples  */
proc optmodel;
   number m = 3, n = 4;
   var x{1..4} init 1;
   string y = 'c';
   put ({<'a', x[3]>, <'b', m>, <y, m/n>});
quit;

/*  SETOF Aggregation Expression  */
proc optmodel;
   put (setof{i in 1..3}<i, i*i, i**3>);
quit;

/*  SLICE Expression  */
proc optmodel;
   put (slice(<1,*>, {<1,3>, <1,0>, <3,1>}));
   put (slice(<*,2,*>, {<1,2,3>, <2,4,3>, <2,2,5>}));
quit;

/*  Warshall's Algorithm  */
proc optmodel;
   set<str,str> dep = {<'B','A'>, <'C','B'>, <'D','C'>};
   set<str,str> cl;
   set<str> cn;
   cl = dep;
   cn =  (setof{<i,j> in dep} i) inter (setof{<i,j> in dep} j);
   for {node in cn}
       cl = cl union (slice(<*,node>,cl) cross slice(<node,*>,cl));
   put cl;
quit;

/*  SUM Expression   */
proc optmodel;
   put (sum {i in 1..10} i);  /* outputs 55 */
quit;

/*  SYMDIFF Expression  */
proc optmodel;
   put ({1,3} symdiff {2,3}); /* outputs {1,2} */
quit;

/*  Tuple Expression  */
proc optmodel;
   put (<1,2,3> in setof{i in 1..2}<i,i+1,i+2>);
   put ({<1,'a'>, <2,'b'>} cross {<3,'c'>, <4,'d'>});
quit;

/*  Tuple Comparison  */
proc optmodel;
   put (<1,'a'> = <1,'a'>)=;  /* true, all elements equal */
   put (<1,'a'> = <1,'b'>)=;  /* false, at least one element unequal */
   put (<1,'a'> ~= <1,'b'>)=; /* true, at least one element unequal */
   put (<1,2,5> < <1,3,4>)=;  /* true, compared left to right, 2 < 3 */
quit;

/*  UNION Expression  */
proc optmodel;
   put ({1,3} union {2,3}); /* outputs {1,3,2} */
quit;

/*  UNION Aggregation Expression  */
proc optmodel;
   put (union{i in 1..3} i..i+3); /* outputs {1,2,3,4,5,6} */
quit;

/*  WITHIN Expression  */
proc optmodel;
   put ({1,3} within {2,3});     /* outputs 0 */
   put ({1,3} not within {2,3}); /* outputs 1 */
   put ({1,3} within {1,2,3});   /* outputs 1 */
quit;

/*  Data Set Input/Output  */
proc optmodel;
   set LOCS = {'New York', 'Washington', 'Boston'};  /* locations */
   set DOW = 1..7;  /* day of week */
   var s{LOCS, DOW} init 1;
   create data soldata from [location day_of_week]={LOCS, DOW} sale=s;
quit;

data pdat;
   input p $ maxprod cost;
   datalines;
ABQ    12  0.7
MIA     9  0.6
CHI    14  0.5
;

proc optmodel;
   set<string> plants;
   var prod{plants} >= 0;
   number cost{plants};
   read data pdat into plants=[p] prod.ub=maxprod cost;
quit;

data abcd;
   input letter $ @@;
   datalines;
a b c d
;

/*  READ DATA Statement: Explicit Indexing  */
proc optmodel;
   set<num> subscripts=1..4;
   string  letter{subscripts};
   read data abcd into [_N_] letter[5-_N_];
   print letter;
quit;

data revdata;
   input month rev @@;
   datalines;
1 200 2 345 3 362 4 958
5 659 6 804 7 487 8 146
9 683 10 732 11 652 12 469
;

/*  CREATE DATA Statement: Explicit Indexing  */
proc optmodel;
   set m = 1..3;
   var revenue{1..12};
   read data revdata into [_N_] revenue=rev;
   create data qtr1 from [month]=m revenue[month];
   create data qtr2 from [month]=m revenue[month+3];
   create data qtr3 from [month]=m revenue[month+6];
   create data qtr4 from [month]=m revenue[month+9];

/*  CREATE DATA Statement: Dynamic Data Set Names  */
proc optmodel;
   set m = 1..3;
   var revenue{1..12};
   read data revdata into [_N_] revenue=rev;
   for {q in 1..4}
      create data ("qtr" || q)
             from [month]=m revenue[month+(q-1)*3];

/*  Formatted Output  */
proc optmodel;
   number a=1.7, b=2.8;
   set s={a,b};
   put a b;        /* list output */
   put a= b=;      /* named output */
   put 'Value A: ' a 8.1 @30 'Value B: ' b 8.; /* formatted */
   string str='Ratio (A/B) is:';
   put str (a/b);  /* strings and expressions */
   put s=;         /* named set output */
quit;

/*  PRINT Statement: List Form Output  */
proc optmodel;
   number square{i in 0..5} = i*i;
   number recip{i in 1..5} = 1/i;
   print square recip;

/*  PRINT Statement: Matrix Form Output  */
proc optmodel;
   set R=1..6;
   set C=1..4;
   number a{i in R, j in C} = 10*i+j;
   print a;
quit;

/*  ODS Table and Variable Names  */
proc optmodel printlevel=2;
   ods output PrintTable=expt ProblemSummary=exps DerivMethods=exdm
              SolverOptions=exso SolutionSummary=exss OptStatistics=exos
              Timing=exti;
   var x{1..2} >= 0;
   min z = 2*x[1] + 3 * x[2] + x[1]**2 + 10*x[2]**2
           + 2.5*x[1]*x[2] + x[1]**3;
   con c1: x[1] - x[2] <= 1;
   con c2: x[1] + 2*x[2] >= 100;
   solve;
   print x;
quit;

title2 'PrintTable';
proc print data=expt;
run;

title2 'ProblemSummary';
proc print data=exps;
run;

title2 'SolverOptions';
proc print data=exso;
run;

title2 'DerivMethods';
proc print data=exdm;
run;

title2 'SolutionSummary';
proc print data=exss;
run;

title2 'OptStatistics';
proc print data=exos;
run;

title2 'Timing';
proc print data=exti;
run;

title2;

/*  Constraints  */
proc optmodel;
   var x, y;
   min r = x**2 + y**2;
   con c: x+y >= 1;
   solve;
   print x y;
quit;

/*  Expanding a Standardized Constraint  */
proc optmodel;
   var x{1..3};
   con b: sum{i in 1..3}(x[i] - i) = 0;
   expand b;
quit;

/*  Expansion and Body Values for a Standardized Constraint  */
proc optmodel;
    var x init 1;
    con c1: x**2 <= 5;
    con c2: 5 >= x**2;
    con c3: -x**2 >= -5;
    con c4: -5 <= -x**2;
    expand;
    print c1.body c2.body c3.body c4.body;
quit;

/*  Constraints: Dual Values  */
proc optmodel;
   var x;
   min o1 = x**2;
   con c1: 1 <= x <= 2;
   solve;
   print (c1.dual);

   max o2 = -x**2;
   solve;
   print (c1.dual);
quit;

/*  Suffixes  */
proc optmodel;
   var v >= 0 <= 2 init 1 suffixes=(label='X0000000' status init 'L');

   num n init 10;
   num ubound{1..N} init 2;
   var x{i in 1..N} >= 0 <= ubound[i] suffixes(label='X'||put(i,z7.));
   expand;
   print _var_.name _var_.label;
quit;

proc optmodel;
   var x, y;
   min z = x + y;
   con c: x + 2*y <= 3;
   solve;
   print x.lb x.ub x.status x.sol;
   print y.lb y.ub y.status y.sol;
   print c.lb c.ub c.body c.dual;

   x.lb=0;
   y.lb=0;
   c.lb=1;
   solve;
   print x.lb x.ub x.status x.sol;
   print y.lb y.ub y.status y.sol;
   print c.lb c.ub c.body c.dual;
quit;

/*  Dual Value  */
proc optmodel;
   var x, y;
   min z = 6*x + 7*y;
   con
      4*x +   y >=  5,
       -x - 3*y <= -4,
        x +   y <=  4;
   solve;
   print x y;
   expand _ACON_ ;
   print _ACON_.dual _ACON_.body;

   _ACON_[1].lb = _ACON_[1].lb - 1;
   solve;
   _ACON_[2].ub = _ACON_[2].ub + 1;
   solve;
quit;

/*  Reduced Costs  */
proc optmodel;
   var x >= 0, y >= 0, z >= 0;
   max cost = 4*x + 3*y - 5*z;
   con
      -x + y + 5*z <= 15,
      3*x - 2*y - z <= 12,
      2*x + 4*y + 2*z <= 16;
   solve;
   print x y z;
   print x.rc y.rc z.rc;
quit;

/*   Presolver  */
proc optmodel;
   var x >= 3;
   var y >= 0;
   con c: x + y <= 7;
   con c2: x + 2*y <= 15; /* redundant */
   expand/solve;
quit;

/*  Model Update  */
proc optmodel;
   var x, y;
   min r = x**2 + y**2;
   con c: x+y >= 1;
   solve;
   print x y;

   drop c;
   solve;
   print x y;

   restore c;
   fix x=0.3;
   solve;
   print x y c.dual;

   unfix x;
   solve;
   print x y c.dual;
quit;

/*  Multiple Solutions  */
proc optmodel printlevel=0;
   var x{1..2} integer >= 1 <= 3;
   con c: alldiff(x);
   solve with clp / allsolns;
   print _NSOL_;
   print {j in 1..2, i in 1.._NSOL_} x[j].sol[i];
   create data solout from [sol]={i in 1.._NSOL_}
       {j in 1..2} <col("x"||j)=x[j].sol[i]> ;

/*  Problem Symbol Output  */
proc optmodel printlevel=0;
   var x1 >= 0, x2 >= 0, x3 >= 0, x4 >= 0, x5 >= 0;

   minimize z = x1 + x2 + x3 + x4;

   con a1: x1 + x2 + x3          <= 4;
   con a2:               x4 + x5 <= 6;
   con a3: x1 +          x4      >= 5;
   con a4:      x2 +          x5 >= 2;
   con a5:           x3          >= 3;

   solve with lp;

   print _var_.name _var_ _var_.rc _var_.status;
   print _con_.name _con_.lb _con_.body _con_.ub _con_.dual _con_.status;

/*  OPTMODEL Options  */
proc optmodel;
   number n = 1/7;
   for {i in 1..9 by 4}
   do;
      reset options pdigits=(i);
      print n;
   end;
   reset options pdigits; /* reset to default */
quit;

/*  Calling FCMP Functions  */
proc fcmp outlib=work.funcs.test;
   /* arithmetic geometric mean */
   function agm(a0, b0);
      a = a0; b = b0;
      if a<=0 or b<=0 then return(0);
      do until( a - b < a/1e12 );
         a1 = 0.5*a + 0.5*b;
         b1 = sqrt(a*b);
         a = a1; b = b1;
      end;
      return(a);
   endsub;
run;

/* libraries must be specified with the CMPLIB option */
option cmplib=work.funcs;

proc optmodel;
   title 'AGM(1,2)';
   print (agm(1,2));

proc optmodel;
   title 'AGM optimization';
   /* find x where agm(1,x) == 23 */
   var x init 1;
   num c = 23;
   min z = (agm(1,x)-c)^2;
   solve;
   print x;
quit;

/*  FCMP Output Arguments  */
proc fcmp outlib=work.funcs.test;
   subroutine do_sqr(x, sq, text $);
      outargs sq, text;
      sq = x*x;
      text = 'This is an example of output arguments';
   endsub;
run;

option cmplib=work.funcs;

proc optmodel;
   string s init repeat(' ', 79); /* reserve 80 bytes */
   number n;
   call do_sqr(7, n, s);
   print s n;
quit;

/*  Constant Matrix Argument to an FCMP Function  */
proc fcmp outlib=work.funcs.test;
   function evalpoly(x, coeff[*]);
      z = 0;
      do i = dim1(coeff) to 1 by -1;
         z = z * x + coeff[i];
      end;
      return (z);
   endsub;
run;

option cmplib=work.funcs;

proc optmodel;
   num coeff{1..3} = [1, -2, 1];
   var x;
   min z=evalpoly(x, coeff);
   solve;
   print x;
quit;

/*  OPTMODEL Declarations for FCMP Matrices  */
proc fcmp outlib=work.funcs.test;
   subroutine mattest(x[*]);
      put x[1]=;
   endsub;
   subroutine mattest2(x[*,*]);
      put x[1,1]=;
   endsub;
run;

option cmplib=work.funcs;

proc optmodel;
   /* the following arrays can be used as matrices */
   num N init 3;
   num mat1{1..N} init 0;
   call mattest(mat1);     /* OK */
   set S1 = 1..5;
   num mat2{S1} init 0;
   call mattest(mat2);     /* OK */
   set S2 = {S1,S1};
   num mat3{S2} init 0;
   call mattest2(mat3);    /* OK */
   num mat4{S1 cross S1} init 0;
   call mattest2(mat4);    /* OK */
   num L init 1;
   num mat5{L..N} init 0;
   call mattest(mat5);     /* OK */
   set S3 init S1;
   num mat6{S3} init 0;
   call mattest(mat6);     /* OK */

   /* some errors are detected at execution time */
   S3 = 2..5;
   call mattest(mat6);     /* ERROR: lower bound not 1 */
   S3 = {1, 3, 4, 5};
   call mattest(mat6);     /* ERROR: missing index 2 */
   L = 0;
   call mattest(mat5);     /* ERROR: lower bound not 1 */

   /* the following arrays cannot be used as matrices */
   num arr1{1..10 by 3};   /* step size is not 1 */
   call mattest(arr1);     /* ERROR */
   num arr2{i in 1..N, j in 1..N: j >= i}; /* selection expression used */
   call mattest2(arr2);    /* ERROR */
   num arr3{i in 1..N, j in 1..i};         /* index dependency on 'i' */
   call mattest2(arr3);    /* ERROR */
quit;

/*  More on Index Sets  */
proc optmodel;
   put (setof {i in 1..3, j in 1..3 : j >= i} <i, j>);
quit;

proc optmodel;
   put ({i in 1..3, i..3});
quit;

/*  SLICE Operator  */
proc optmodel;
   number N = 3;
   set<num,str> S = {<1,'a'>,<2,'b'>,<3,'a'>,<4,'b'>};
   put ({i in 1..N, <(i),j> in S});
   put ({i in 1..N, j in slice(<i,*>, S)});
quit;

/*  Solver Status Parameters  */
proc optmodel printlevel=0;
   var x init 1 >= 0.001;
   min z=sin(x)/x;
   solve;
   for {k in /STATUS SOLUTION_STATUS ALGORITHM/}
      put _OROPTMODEL_STR_[k]=;
   for {k in /OBJECTIVE ITERATIONS/}
      put _OROPTMODEL_NUM_[k]=;
quit;

/*  Macro Resolution  */
proc optmodel;
   var x, y;
   min z=x**2 + (x*y-1)**2;
   for {n in 1..3} do;
      fix x=n;
      solve;
      %put Line 1 &_OROPTMODEL_;
      put 'Line 2 ' (symget("_OROPTMODEL_"));
   end;
quit;

/*  PROC NLP Example 7  */
proc nlp all;
   parms amountx amounty amounta amountb amountc
         pooltox pooltoy ctox ctoy pools = 1;
   bounds 0 <= amountx amounty amounta amountb amountc,
               amountx <= 100,
               amounty <= 200,
          0 <= pooltox pooltoy ctox ctoy,
          1 <= pools <= 3;
   lincon amounta + amountb = pooltox + pooltoy,
          pooltox + ctox = amountx,
          pooltoy + ctoy = amounty,
          ctox + ctoy    = amountc;
   nlincon nlc1-nlc2 >= 0.,
           nlc3 = 0.;
   max f;
   costa = 6; costb = 16; costc = 10;
   costx = 9; costy = 15;
   f = costx * amountx + costy * amounty
       - costa * amounta - costb * amountb - costc * amountc;
   nlc1 = 2.5 * amountx - pools * pooltox - 2. * ctox;
   nlc2 = 1.5 * amounty - pools * pooltoy - 2. * ctoy;
   nlc3 = 3 * amounta + amountb - pools * (amounta + amountb);
run;

proc optmodel;
   var amountx init 1 >= 0 <= 100,
       amounty init 1 >= 0 <= 200;
   var amounta init 1 >= 0,
       amountb init 1 >= 0,
       amountc init 1 >= 0;
   var pooltox init 1 >= 0,
       pooltoy init 1 >= 0;
   var ctox init 1 >= 0,
       ctoy init 1 >= 0;
   var pools init 1 >=1 <= 3;
   con amounta + amountb = pooltox + pooltoy,
       pooltox + ctox = amountx,
       pooltoy + ctoy = amounty,
       ctox + ctoy    = amountc;
   number costa, costb, costc, costx, costy;
   costa = 6; costb = 16; costc = 10;
   costx = 9; costy = 15;
   max f = costx * amountx + costy * amounty
          - costa * amounta - costb * amountb - costc * amountc;
   con nlc1: 2.5 * amountx - pools * pooltox - 2. * ctox >= 0,
       nlc2: 1.5 * amounty - pools * pooltoy - 2. * ctoy >= 0,
       nlc3: 3 * amounta + amountb - pools * (amounta + amountb)
              = 0;
   solve;
   print amountx amounty amounta amountb amountc
         pooltox pooltoy ctox ctoy pools;