ERROR: Unexpected exception while loading model

[Kiran Sheth]

Hi,

I get "ERROR: Unexpected exception while loading model" message while running the attached model in pyomo, However, I am not sure how can I diagnose what is causing the error. I have tried various command line switches (incl. --debug), but have not been able to make any progress. Can someone help me either diagnose the error, or suggest a way to diagnose the error? 

Many "Thanks" in advance.

-Kiran

[Gabriel Hackebeil]

A few things:

(1) You don’t need to use the SolverFactory in a model file that you intend to use with the pyomo command. This is done automatically via the --solver=‘ipopt’ command-line option.

(2) You are not loading external data to construct your instance, so you should be starting from a ConcreteModel object (rather than AbstractModel). If you stick with an AbstractModel, all of your Constraint declarations need to be switched over to using rules like your Objective declaration.

(3) You basically already have a script that you can run without the pyomo command (e.g., python LP.py, using the coopr-installed python executable). You don’t need a pyomo_postprocess declaration to print your results in this case. You can just print the name and value of whatever components you want to inspect on your model after calling instance.load(results):

print instance.CW.cname(), value(instance.CW)

instance.pprint()

display(instance)

print results

(4) Pyomo provides the Block() modeling component for name-spacing other components, which is what I assume you are trying to do with the Struct class. Blocks can be initialized using a rule or modified in a concrete fashion just like a model. E.g.:

model = ConcreteModel()

model.Feed = Block()

model.Feed.Q = Var(…)

Gabe

On Aug 23,

 2014, at 1:11 PM, Kiran Sheth <krshe...@gmail.com> wrote:

<LP.py>

[Kiran Sheth]

Gabriel,

Many "Thanks" for your assistance.

(1) I am aware of the --solver switch. I was just being lazy to not type it at the command line, so decided to put it in the script instead.

(2) I am a novice user of Pyomo, and not clear on the difference between Abstract Vs Concrete model. But, my ultimate model will have something like

        instance1 = model.create()

        instance1.Feed.Q = ...

        ...

        instance1.preprocess()

        instance2 = model.create()

        instance2.Feed.Q = ..

        ...

        instance2.preprocess()

where instance1 and instance2 will have same model structure, but different data. So, my naive thinking was that I should declare AbstractModel with correct structure once, and invoke it as many times as I want using model.create(). Any guidance you can offer to accomplish the goal, is appreciated.

(3) I was switching back and forth between pyomo and python to try to debug the problem I was having. So, I left the pyomo_postprocess method in the script just being lazy.

(4) Using Block() instead of the Struct() class certainly allows both pyomo or python to load and attempt to solve the model. So, now I have few related questions.

Q(4.1) I downloaded pyomo.pdf and CooprGettingStarted.pdf from the web. Neither one mentions anything about Block component. Is there a better place to learn all the features of pyomo?

Q(4.2) I was using Stuct class to declare a group of variables (Q,T,Ca,Cb,Cc), and allowing model.Feed and model.Prod to inherit the group. Using block component, now I have to define

        model.Feed.Q = ...

        model.Feed.T = ...

        ...

        model.Prod.Q = ...

        model.Prod.T = ...

        ...

etc.   i.e. all the variables need to be declared twice. With Struct I just had to declare it once. Am I missing something?

Q(4.3) I got the definition of Struct class from the web. Obviously since it is not working, something is wrong with it. By any chance, do you see something that sticks out as wrong?

Q(4.4) With Block() now the model loads, and as mentioned earlier, attempts to solve. But it seems that it is not solved using the values I provided via instance.Feed.Q = ... commands. I thought that once I assign values via these commands, and issue instance.preprocess command, the model will solve using these values. Is there any other command that I need to issue?

The modified script file is attached. Again, appreciate your assistance.

-Kiran

[Gabriel Hackebeil]

(2) I am a novice user of Pyomo, and not clear on the difference between Abstract Vs Concrete model. But, my ultimate model will have something like

        instance1 = model.create()

        instance1.Feed.Q = ...

        ...

        instance1.preprocess()

        instance2 = model.create()

        instance2.Feed.Q = ..

        ...

        instance2.preprocess()

where instance1 and instance2 will have same model structure, but different data. So, my naive thinking was that I should declare AbstractModel with correct structure once, and invoke it as many times as I want using model.create(). Any guidance you can offer to accomplish the goal, is appreciated.

Calling model.create() is only necessary if you are using an AbstractModel. Note that if you are starting from a ConcreteModel, model.create() is a no-op and returns a reference to itself. If you want a copy of a ConcreteModel, you can call model.clone(). This post (https://groups.google.com/forum/?hl=en#!topic/coopr-forum/kVfC1pZmhqo) gives a nice explanation of the differences between these two modeling paradigms. In either case, it is still appropriate to call preprocess() after any major changes like adding constraints, fixing variables, updating parameter values, etc.

Q(4.1) I downloaded pyomo.pdf and CooprGettingStarted.pdf from the web. Neither one mentions anything about Block component. Is there a better place to learn all the features of pyomo?

We’re past due on adding documentation for this component. I’ll add a ticket for that.

Q(4.2) I was using Stuct class to declare a group of variables (Q,T,Ca,Cb,Cc), and allowing model.Feed and model.Prod to inherit the group. Using block component, now I have to define

        model.Feed.Q = ...

        model.Feed.T = ...

        ...

        model.Prod.Q = ...

        model.Prod.T = ...

        ...

etc.   i.e. all the variables need to be declared twice. With Struct I just had to declare it once. Am I missing something?

As with other components, you can initialize a Block using a rule. E.g.,

def stream_rule(block):

   block.Q = Var(…)

   ...

   block.CC = Var(…)

model.Feed = Block(rule=stream_rule)

model.Prod = Block(rule=stream_rule)

Alternatively, you can define a function that creates a block for you:

def myblock(model):

   b = Block()

   b.x = Var()

   # Assuming model is a ConcreteModel, otherwise

   # a rule would be required here

   b.c = Constraint(expr=b.x == model.v)

   return b

model.b = myblock(model)

Q(4.3) I got the definition of Struct class from the web. Obviously since it is not working, something is wrong with it. By any chance, do you see something that sticks out as wrong?

A Pyomo Block maintains a specialized list of the Pyomo components that have been assigned to it, including other Blocks (AbstractModel and ConcreteModel are subclasses of Block). This way, when a solver interface wants to collect the set of constraints on a model, it can do so using code similar to the following:

for block in model.all_blocks():
   for var in block.active_subcomponents(Constraint).values()

      ...

A Pyomo model wouldn’t recognize you’re Struct class as being a Block, so it obviously would be ignored by the above loop.

Q(4.4) With Block() now the model loads, and as mentioned earlier, attempts to solve. But it seems that it is not solved using the values I provided via instance.Feed.Q = ... commands. I thought that once I assign values via these commands, and issue instance.preprocess command, the model will solve using these values. Is there any other command that I need to issue?

In the model you provided, you have Feed.Q defined as a variable. Therefore, a statement like "instance.Feed.Q =  0.25” is actually updating the current value of the variable on your model, which in turn gets used as the starting point for a solver like Ipopt. This can indeed result in a different solution being returned for nonlinear models that are not convex, but I assume you are really expecting Feed.Q to be fixed to the value you are assigning to it. In this case you should use something like “instance.Feed.Q.fix(0.25)”.

If you’re are always expecting this kind of behavior it might be more appropriate to redefine some of these components as parameters rather than variables. E.g.,

model.Feed.Q = Param(mutable=True, initialize=0.3)

…

model.preprocess()

… solve

… load results

model.Feed.Q = 0.25

model.preprocess()

… solve

… load results

Using the “mutable=True” keyword with Param allows its value to be updated in any expressions where it is used after the model in constructed. The default behavior for a Param is to be immutable and so will raise an exception if you try update its value after construction. The online documentation is a bit sparse for this feature as well.

Gabe

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

[Kiran Sheth]

Gabriel,

This has been of great help. I have been able to make good progress,

Now I have a ConcreteModel LP_RXN, and another ConcreteModel LP_MXR. Both these models have their own constraints that are written as 

    LP_RXN.CNSTRx = Constraint(expr=(...)) and 

    LP_MXR.CNSTRx = Constraint(expr=(...)) respectively.

I have also created instances of both LP_RXN (e.g. RX) and LP_MXR (e.g. MX), and have been able to successfully run them with appropriate data.

However, I need to build, what we in process industry call a flowsheet, where the elements of Prod block of RX are "connected" to the Feed block of MX. I conceptually want to write additional constraints such as

Cxn_Q = Constraint(expr=(0 == MX.Feed.Q - RX.Prod.Q))

However, since this constraint does not belong to either of the two RX or MX models, I am not sure how to write it. Is there a way to write such global / orphan constraints in pyomo?

Also, since RX.Prod and MX.Feed blocks have multiple elements (Q, T, Ca, Cb, Cc), is there a clever way to loop over all the elements to connect all the elements?

[Gabriel Hackebeil]

However, since this constraint does not belong to either of the two RX or MX models, I am not sure how to write it. Is there a way to write such global / orphan constraints in pyomo?

Every Pyomo component should live on exactly 1 model/block. If you want to create a new “master” model out of other independent models, the workflow would look something like:

LP_RXN = ConcreteModel()

… add components

… solve / load results

LP_MXR = ConcreteModel()

… add components

… solve / load results

# Free any variables that were fixed in order to warmstart the sub model solutions

LP_RXN.Q.free()

LP_MXR.Q.free()

...

# Create the master flowsheet model

master = ConcreteModel()

master.RXN = LP_RXN

master.MXR = LP_MXR

# Add master model objective

master.OBJ = Objective(…)

# Make sure we only have one active Objective

master.RXN.OBJ.deactivate()

master.MXR.OBJ.deactivate()

# Add linking constraints

master.Cxn_Q = Constraint(expr=master.RXN.Q == master.MXR.Q)

…

master.preprocess()

… solve / load results

Note that once you add a model onto another master model, it cannot be solved independently until it is removed from the master. E.g.,

master.del_component(‘RXN’)

LP_RXN.OBJ.activate()

LP_RXN.preprocess()

… solve

You should also get some kind of warning if you add a sub model to a different master model before removing it from the first (via the del_component() method).

Also, since RX.Prod and MX.Feed blocks have multiple elements (Q, T, Ca, Cb, Cc), is there a clever way to loop over all the elements to connect all the elements?

We have implemented a “Connector” component which was designed for this kind of modeling situation. However, I believe this is still a bit in the experimental phase (someone might correct me here). I can’t find any online documentation, but I do see an example in coopr.pyomo/examples/pyomo/connectors/ that may be of use.

I would, however, recommend first just defining a function for “connecting” different process units before trying to experiment with the Connector component (to avoid any potential debugging headaches). E.g.,

def connect(unit1, unit2):

b = Block()

b.Q_link = Constraint(expr=unit1.Q == unit2.Q)

b.T_link = Constraint(expr=unit1.T == unit2.T)

…

return b

master.RXN_MXR_connect = connect(master.RXN, master.MXR)

If you do try out the Connector component, any suggestions for improving / expanding the interface is welcome.

Gabe

[Kiran Sheth]

Gabriel,

I don't know how to thank you. But your advise was spot on. While I have not tried the connectors, using the connect method below, I have been able to create the flowsheet (master) with several units (RX1, MXR, RX2) connected to each other, and optimize the flowsheet.

Now I am trying to generate a report. I had this fragment of code (shamelessly stolen from pyomo / coopr documentation) working at individual model level (i.e. when there was no flowsheet / master).

def myreport (options, instance, results):

    from coopr.pyomo import Var

    for v in instance.active_components(Var):

        varobject = getattr(instance, v)

        for index in varobject:

            print "    ",index, varobject[index], \

                varobject[index].value, \

                "SV(LB) ", instance.ipopt_zL_out.getValue(varobject[index]), \

                "SV(UB) ", instance.ipopt_zU_out.getValue(varobject[index])

        #

    #

However, it does not work at the flowsheet (master) level. Conceptually, I was printing shadow values of all variables using RX1.ipopt_z{L|U}_out.getValue(RX1.Prod.T) type command. But, now I need to issue LP.ipopt.z{L|U}_out.getValue(LP.RX1.Prod.T) command. FYI - LP is the name I used for master. I suspect that I need two loops - outer one for accessing underlying units (RX1, MXR, RX2) and then inner one for the variables in that unit. However, I can not figure out how to get the list of underlying units in the master. 

Furthermore, I need to print optimization variables that the solver optimized. I designated them as free in my model. In that case, free variables include optimization variables that were used for optimizing the objective function, as well as dependent variables that were needed for satisfying the model constraints. Is there another variable attribute I should be using for differentiating optimization variables?

-Kiran

[Gabriel Hackebeil]

However, it does not work at the flowsheet (master) level. Conceptually, I was printing shadow values of all variables using RX1.ipopt_z{L|U}_out.getValue(RX1.Prod.T) type command. But, now I need to issue LP.ipopt.z{L|U}_out.getValue(LP.RX1.Prod.T) command. FYI - LP is the name I used for master. I suspect that I need two loops - outer one for accessing underlying units (RX1, MXR, RX2) and then inner one for the variables in that unit. However, I can not figure out how to get the list of underlying units in the master. 

There are many ways to loop over subblocks on a model (or Block) and there’s subtle differences to be aware of when you have indexed vs. non-indexed Blocks as well as multi-level nesting of sub-blocks. I’ve attached a small example script that should help you wrap your head around things. FWIW, I think you are probably looking for something like “LOOP 3” in the example script.

Another simple approach might be to create an indexed Block that refers to the process units you are connecting. E.g.,

RXN = ConcreteModel()

model.UnitNames = Set()
model.Units = Block(model.UnitNames)

model.UnitNames.add('RXN')
model.Units['RXN'].m = RXN

for unit_name in model.UnitNames:
    unit_model = model.Units[unit_name].m
    print unit_model.cname(True)

Furthermore, I need to print optimization variables that the solver optimized. I designated them as free in my model. In that case, free variables include optimization variables that were used for optimizing the objective function, as well as dependent variables that were needed for satisfying the model constraints. Is there another variable attribute I should be using for differentiating optimization variables?

I don’t think there is really a way for Pyomo to distinguish between variables this way.  Although you could partition the variables on your unit models by placing them on separate subblocks, that way you could generically loop over the “optimization” variables designated for a given unit by looking in the “optimization” block. Although, creating a list of variables or variable names and placing that on the unit model might be just as effective.

BTW, this thread is turning into a good template for some initial online Block documentation, so don’t hesitate to ask if you’re unclear on anything Block related.

Gabe

[Kiran Sheth]

Gabriel,

Looping over LP.all_blocks() helped a lot. Now, I have been able to print the report that I had in mind. Couple of related questions-

(1) I was using insrance.active_components(Var), instance.active_components(Param) and instance.active_components(Constraint) commands, and your script also has instance.active_components(Block) command. Which begs the question, what other arguments the instance.active_components command takes? 

(2) Right now all my constraints are written as 

    Inst.Cnstr1 = Constraint(expr = ...)

    Inst.Cnstr2 = Constraint(expr = ...)

I am also familiar with

    def constraint_rule(model, i):

        # return the expression for the constraint for i

        return sum(model.a[i,j] * model.x[j] for j in model.J) >= model.b[i]

    model.LIEQ = Constraint(model.I, rule=constraint_rule)

form of writing a set of constraints by looping over a RangeSet. 

However. I am not sure if I can write something like

    def all_constraint_rule(model):

        return (Cnstr1 = ... , Cnstr2= ... , ....)

    model.LIEQ = Constraint(model, rule=all_constraint_rule)

or

    def all_constraint_rule(model):

        Cnstr1 = ...

        Cnstr2 = ....

        return (Cnstr1, Cnstr2, ....)

    model.LIEQ = Constraint(model, rule=all_constraint_rule)

so that I do not have to write all my constraints as expressions, and I can use the constraint_rule instead.

-Kiran

[Gabriel Hackebeil]

Again, sorry for the late response.

(1) I was using insrance.active_components(Var), instance.active_components(Param) and instance.active_components(Constraint) commands, and your script also has instance.active_components(Block) command. Which begs the question, what other arguments the instance.active_components command takes? 

Set, Param, Var, Objective, Block, Constraint, Expression, as well as a few other less common “Component” types found in Coopr.

(2) Right now all my constraints are written as 

    Inst.Cnstr1 = Constraint(expr = ...)

    Inst.Cnstr2 = Constraint(expr = …)

    ...

You are probably looking for something like the ConstraintList object. This is just a Constraint that implicitly keeps track of a dynamic index (meaning it appears in the list returned by active_components(Constraint)). You can provide a python generator to initialize it with a set of constraint expressions

def c_rule(model):
    yield model.x <= 1
    yield 1 <= model.y + model.x <= 10
    yield model.y**2 + model.x == 1
model.c = ConstraintList(rule=c_rule)

or just add individual constraints to its collection later on with

model.c.add(model.x == 0)

Gabe