Constraint with If statement, variables and value()

[Thomas Frei]

Hi  (sorry if this gets posted twice, I had a small mistake in my first post) 

I'm trying to use a constraint in which an if statement on a variable is present. It seems that this has to be done with the value() function. 

Unfortunately, I get weird model behaviour, when I do that. After I initialized the variable in order to use value() the result is that pyomo stays

inside the corresponding if clause even when p_diff takes on values for which the if statement isn't true anymore and one of the other conditions should be used.

In the code below for example pyomo keeps on skipping the constraint, since p_diff was initialized with 0. If I initialize it with a non-negative value it stays in 

one of the other if statements until **0.5394 can't be calculated anymore and I get an error

I have tried multiple ways, but nothing seems to work. Does anyone know how to handle such a case? 

Example:

m.p_diff=Var(m.a, m.t, initialize=0, doc="Pressure difference for panhandle equation")

def p_diff_def(m, a, t):
    return m.p_diff[a,t]==(((10**3*m.p[m.arc_n0[a],t])**2-(10**3*m.p[m.arc_n1[a],t])**2)/(m.spec_grav**0.8539*m.T_avg*m.L[a]*m.z_gas))
m.p_diff_Con=Constraint(m.a, m.t, rule=p_diff_def)


###Panhandle A equation
def panhandle_def(m, a, t):
    if value(m.p_diff[a,t])>=0 and value(m.p_diff[a,t])!=0:
                    return m.m_left_p[a,t]*10**3/m.rho_norm*24==4.5965*10**-3*m.e_pipe*(m.T_std/(m.p_atm*10**3))**1.0788\
                                                            *m.p_diff[a,t]**0.5394*m.D[a]**2.6182
    if value(m.p_diff[a,t])<=0 and value(m.p_diff[a,t])!=0:
                    return m.m_right_p[a,t]*10**3/m.rho_norm*24==4.5965*10**-3*m.e_pipe*(m.T_std/(m.p_atm*10**3))**1.0788\
                                                            *(-m.p_diff[a,t])**0.5394*m.D[a]**2.6182
    else:
                    return Constraint.Skip
m.panh_con=Constraint(m.a, m.t, rule=panhandle_def, doc="Panhandle A equation")

[David Woodruff]

Here's a short version of an answer:

The *if* statements in constraint construction are evaluated when the model is instantiated. The *if* is in no way, shape, or form passed to the solver. So only the initialized value of the variable matters to this *if*.

End of Pyomo discussion.

In general, if you want to "select" constraints based on variable values during the solve, there are standard tricks for this that involve setting up indicator variables (e.g., named delta) and then using them on the right hand side  in BigM expressions such as (1-delta)*BigM....

See the modeling book by HP Williams or other modeling texts for examples and discussion.

  Dave

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Gabriel Hackebeil]

Adding to Dave’s response...

Generally, situations like this require one to introduce discrete variables into the model to represent logical conditions.

However, this is not always necessary. If you use a solver through the NL file interface (which likely you are), there is support for embedding logical conditions within the expression. In Pyomo, this is accomplished by representing the logical condition with an Expr_if declaration. E.g.,

def c_rule(m, i):

    return Expr_if(IF=(m.x <= 0), THEN=(-m.x**3), ELSE=(m.x**3)) * 5 * m.y**2 + ... >= ...)

I don’t know that this is documented yet, so you might need to add an extra import to use it (from pyomo.core.base.expr import Expr_if).

In most situations, this is the wrong way to go because the solver will not be aware that you are expressing logical conditions. However, it can be a useful way (and the only way) to express certain continuous / smooth nonlinear functions. You need to be careful though. If the characteristics of function you end up defining don’t meet the set of assumptions the solver makes (e.g., twice continuously differentiable), you can end of thinking you have a good solution when the solver just returned garbage.

I recall at one point working on a model that used a more simple pressure drop equation where this was applicable (e.g., p_diff^2). Unfortunately, I don’t think this approach is applicable in your case (unless you do some smoothing around zero). You could imagine doing some polynomial approximation near zero, and embedding a few of the Expr_if statements to define a smooth, differentiable function.

Regards,

Gabe

[Thomas Frei]

Of course!! Thank you Dave and Gabe. This makes sense. As you see can see I'm new to mathematical modelling :) 

I will change all my if clauses to sth like that, with m.p_dir being the flow direction (1=+ or 0=-). I guess then it should work :

###Free Flow direction
def flow_p_direction(m, a, t):
        return m.p[m.arc_n0[a],t]*(m.p_dir[a,t])+(1-m.p_dir[a,t])*m.p[m.arc_n1[a],t]\
               >=m.p[m.arc_n1[a],t]*m.p_dir[a,t]+(1-m.p_dir[a,t])*m.p[m.arc_n0[a],t]+0.0001
m.f_direct_p=Constraint(m.a, m.t, rule=flow_p_direction, doc="Definition of flow direction in passive lines")

###Panhandle A equation
def panhandle_def(m, a, t):

       return (m.p_dir[a,t])*m.m_left_p[a,t]*10**3/m.rho_norm*24+(1-m.p_dir[a,t])*m.m_right_p[a,t]*10**3/m.rho_norm*24\
                  ==4.5965*10**-3*m.e_pipe*(m.T_std/(m.p_atm*10**3))**1.0788\
                    *(((10**3*((m.p_dir[a,t])*m.p[m.arc_n0[a],t]+(1-m.p_dir[a,t])*m.p[m.arc_n1[a],t]))**2
                      -(10**3*((1-m.p_dir[a,t])*m.p[m.arc_n0[a],t]+(m.p_dir[a,t])*m.p[m.arc_n1[a],t]))**2)
                    /(m.spec_grav**0.8539*m.T_avg*m.L[a]*m.z_gas))**0.5394*m.D[a]**2.6182

m.panh_con=Constraint(m.a, m.t, rule=panhandle_def, doc="Panhandle A equation")

[Thomas Frei]

Small follow up question I was hoping you can help me with: I now redefined the if statements. But there is still something weird going on. In the code below I defined that if the flow direction is from left to right (m.p_dir=1) the pressure in the left node should be higher than in the right one. In the second equation then the difference is calculated from the higher pressure minus the smaller one. The strange thing is, that I get an error: "can't evaluate pow' (-57.91, 0.5394) which would mean that the difference turned out to be negative. Do you know why? Is this the case you talked about Gabe where the polynomial expressions have to be used?

###Pressure levels at line ends

def flow_p_direction(m, a, t):
        return m.p[m.arc_n0[a],t]*(m.p_dir[a,t])+(1-m.p_dir[a,t])*m.p[m.arc_n1[a],t]\

                >=m.p[m.arc_n1[a],t]*m.p_dir[a,t]+(1-m.p_dir[a,t])*m.p[m.arc_n0[a],t]+0.001
m.f_direct_p=Constraint(m.a, m.t, rule=flow_p_direction, doc="Definition of pressure at line ends")

###Panhandle A equation
def panhandle_def(m, a, t):
       return (m.p_dir[a,t])*m.m_left_p[a,t]*10**3/m.rho_norm*24+(1-m.p_dir[a,t])*m.m_right_p[a,t]*10**3/m.rho_norm*24\
                  ==4.5965*10**-3*m.e_pipe*(m.T_std/(m.p_atm*10**3))**1.0788\
                    *(((10**3*((m.p_dir[a,t])*m.p[m.arc_n0[a],t]+(1-m.p_dir[a,t])*m.p[m.arc_n1[a],t]))**2
                      -(10**3*((1-m.p_dir[a,t])*m.p[m.arc_n0[a],t]+(m.p_dir[a,t])*m.p[m.arc_n1[a],t]))**2)
                    /(m.spec_grav**0.8539*m.T_avg*m.L[a]*m.z_gas))**0.5394*m.D[a]**2.6182
m.panh_con=Constraint(m.a, m.t, rule=panhandle_def, doc="Panhandle A equation")

[Gabriel Hackebeil]

The strange thing is, that I get an error: "can't evaluate pow' (-57.91, 0.5394) which would mean that the difference turned out to be negative. Do you know why? Is this the case you talked about Gabe where the polynomial expressions have to be used?

This isn’t so strange if you consider that most optimization software only deals with real numbers. The logical conditions you are trying to enforce are not working out the way you expect and at some point the solver steps into negative territory for x^0.5. So this isn’t directly related to what I suggested, but it may still be important to deal with as you are defining a function that is not differentiable at zero.

I’m assuming you are using an MINLP solver and you just haven’t got your logical conditions worked out. It would help to get these working on a very simplified linear version of the model before anything else.

Gabe

[Thomas Frei]

Hi Gabe

Thank you for your answer. I followed your advice and highly simplified the model. Right now basically the only thing that's left is this non-linear pressure relation and I can be sure now that the strange behaviour is not due to a constraint missmatch. The solver is however, only able to find a solution, if I initialize the pressures. Otherwise I get an error: "can't evaluate pow(0, 0.5394)". I was hoping that I could do this without pressure initialization...but it looks like the solver  doesn't consider my second constraint. I tried Couenne, Bonmin and IPOPT, but I guess they all try to do the same. Attached you can find the simplified model. In case you know a different way, please let me know.

Best,

Thomas

[Carl Laird]

Hi Thomas,

If you don’t initialize the pressures yourself, then they will be initialized to zero by default. As you indicated, this is not really a good initialization for the problem you have. When solving nonlinear problems, initialization to some point is generally required.

Regards,

Carl D. Laird
Associate Professor, School of Chemical Engineering, Purdue University
Ph: 765.494.0085

Email: carl...@purdue.edu

<BE_simple.dat><Test_Model.py>

[Thomas Frei]

Hi Carl

Okay, I will do this then. Thank you very much for your help!

Best,

Thomas

[Carl Laird]

No problem.

Regards,

Carl D. Laird
Associate Professor, School of Chemical Engineering, Purdue University
Ph: 765.494.0085

Email: carl...@purdue.edu