topic Re: substitution of variables in Mathematical Optimization, Discrete-Event Simulation, and OR

substitution of variables

emil — Thu, 21 Mar 2013 08:50:38 GMT

Hi everybody.

Look at this simple optmodel statement as an example:

proc optmodel printlevel=0 FD=CENTRAL;

var j;

max s = j;

con c: (if 1<j<=3 then 3 else 5) <=4;

solve;

create data output from j;

quit;

If you execute it, SAS doesn't find a maximum for j (it reaches 1.0001E19 with a default number of iterations), while to me it might be j=3, since the constraint is true only if 1<j<=3. I see that j disappears from the constraint, in the end, so there is no inequality in j, but in any case I don't understand why I cannot write this.

Why?

It looks to me like SAS translates the problem into something with the following 2 constraints:

con c1: i <= 4;

con c2: i = if 1<j<=3 then 3 else 5;

Now, for every j, this is a system made of an equality (i=3 or i=5) and an inequality (i<=4). Since there is no special priority checking the constraints, SAS will hardly find a valid solution: the only two valid solutions are j=3 and j=5. Is this problem equivalent to the first in SAS? Why?

How can I solve it?

Thank you.

Re: substitution of variables

emil — Thu, 21 Mar 2013 09:42:17 GMT

Maybe this is clearer. Here I try to maximize i, but it doesn't give the right solution anyway.

It gives i=4.0000013216, hence 4, while it might be, in my opinion, 3, because of the c1 constraint:

proc optmodel printlevel=0 FD=CENTRAL;

var i,j;

max s = i;

con c1: i = if 1<j<=3 then 3 else 5;

con c2: i<=4;

solve;

create data output from i j;

quit;

Re: substitution of variables

LeoLopes — Thu, 21 Mar 2013 15:12:24 GMT

What you are trying to do with the model? That might help you get a more specific answer. The 1 < j <= 3 snippet indicates to me that you might be trying to do something beyond what the solver was designed to handle.

The main improvement you could try is to add functions of j to the constraint c, instead of using the constants 3 and 5. The constants tell the solver whether the problem is feasible or not, but don't tell it what to try to do to make the solution feasible, or keep it feasible. So the solver keeps trying to find a better solution, because as far as it knows, feasibility isn't getting any better or worse.

There are a couple of other important improvements you could try as well:

Try to avoid strict inequalities, < or >. The theory that underlies the solvers assume that constraints are always specified with <=, >=, or =. This is why PROC OPTMODEL will produce a syntax error if a constraint is defined with a strict inequality. Using a conditional, you can get around that, and it might work, but it might also fail unpredictably. In practice, even when the solver produces the right output, you will get values very very close to what you would have gotten with a <= or >= anyway.
Try to always declare bounds for the variables if you can (in this case it appears that you can).

Here is a snippet that incorporates the guidelines above:

proc optmodel;

var j >= 1 <= 3; /* use bounds if available */

max s = j;

/* tell the solver that if j < 1 then it should be increased */

/* warning: this defines an open set -- use with extreme caution -- avoid if possible */

con c: (if 1 < j then j else j/2 ) >= 1;

solve;

create data output from j;

quit;

Re: substitution of variables

emil — Thu, 21 Mar 2013 15:59:33 GMT

Ok, thanks..maybe I'm using the proc optmodel in the wrong way. You've been very clear, but still I don't understand how I can define i=f(j). In other words, I need i to assume only some fixed values, depending on values of j (continue). Is there a way to do this?

For example:

proc optmodel;

var i,j;

min s = i+j;

/* statement I'd like to be able to define, to have i a variable depentent on j */

" i = if j>=0 then 1 else 5;"

/* constraint on j*/

con c1: j>=0;

solve;

create data output from i j;

quit;

I want the solution to be i+j=1+0=1.

Thank you so much!

Re: substitution of variables

LeoLopes — Thu, 21 Mar 2013 16:46:46 GMT

Thanks. I am glad I was helpful.

There are several ways to do represent the relationship between two variables. Which one is the best depends on the size of the problem you are trying to solve and on the kind of relationship you need to express. Here are some rules of thumb:

If the relationship is a lookup table, i.e., there is no simple way to write i as a function of j, and especially if you have many strict inequalities, you might want to look at PROC CLP. From the example in your response, I think this may be worth a look. See if this example is helpful: SAS/OR(R) 12.1 User's Guide: Constraint Programming.

If the relationship contains thresholds or represents fixed costs, you might be able to use integer variables with the MILP solver in OPTMODEL. If your other constraints are linear, then this is the option most likely to scale to very large problems. See if this example is helpful: SAS/OR(R) 12.1 User's Guide: Mathematical Programming

If the relationship is a piecewise function, i.e., it changes at certain points but is continuous and smooth otherwise, then: if the problem is linear, use MILP. If it is nonlinear, you can try the NLP solver, especially with the multistart option.

Good luck!

Re: substitution of variables

emil — Thu, 21 Mar 2013 17:36:09 GMT

Thank you Leo.

I was looking at proc clp too, anyway, in those examples I don't see relationships between a continuos variable (j in my example) and a variable which I want to discretize (i in my example). If you want, I need to define a piecewise function from j in real numbers to i in integers. A priori I don't know j, I can only bound it.

This discretization must be runtime, j varies (and the objective function is in terms of j) but consequently i (discretization of j) varies itself and I have a constraint on it.

In the end, are you telling me that there's no easy way to re-write my last example defining i in terms of j?

Again, thank you so much for your time.

Re: substitution of variables

LeoLopes — Thu, 21 Mar 2013 19:58:12 GMT

No problem, I am glad to help.

[Edit: apologies -- in an effort to simplify the presentation I had a version of the code below that was incorrect.]

There are several ways to express the rule you propose, but which one will produce the best results depends on context. The rule is non-linear and discontinuous. In theory, that prevents the algorithms from doing clever things to produce results faster. But in practice, the algorithms implement many tricks that often produce great results. So you might have to experiment a little to see which algorithm and formulation will perform best in practice.

For example: here is how you might do it if you have only linear expressions:

proc optmodel;

num jRange = 5; /* suppose J can vary between -5 and 5 */

num delta = .1; /* meaningful deviation from 0 -- Don't make too small.*/

var JPosPart >= 0 <= jRange, JNegPart >= 0 <= jRange ;

impvar J = JPosPart - JNegPart;

var IsJNeg binary;

impvar I = 1 + 4*IsJNeg;

/* The bound is >= 1 only if J <= delta */

con ComputeIsJNeg1:

IsJNeg <= (1/delta) * jNegPart ;

/* The right hand is > 0 only if jNegPart is < delta, but is always <= 1 */

con ComputeIsJNeg2:

IsJNeg >= 1 / jRange * ( jNegPart - delta );

/* The right hand is 1 only if jPosPart is 0, but always >= 0 */

con ComputeIsJNeg3:

IsJNeg <= 1 - (1/(jRange + delta)) * jPosPart;

max s = 10*I + J ;

solve ;

print I J IsJNeg JNegPart JPosPart;

quit;

There are simpler ways if you make assumptions about objective sense. You might just have to experiment with some examples.