Pyomo Homework 3

Pyomo Homework 3#

# Import the libraries you need here for the assignment

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, '../')
    import helper
helper.set_plotting_style()

# Add your solution here

Pyomo.DAE: Reaction Kinetics#

Consider the chemical reaction $$A \Leftrightarrow B \Leftrightarrow C$$

which is modeling with the following differential algebraic equations:

\[\begin{split}\begin{align*} \frac{dz_A}{dt} &= -p_1 z_A(t) + p_2 z_B(t),\quad z_A(0)=1 \\ \frac{dz_B}{dt} &= p_1 z_A(t) - (p_2 + p_3) z_B(t) + p_4 z_C(t), \quad z_B(0)=0 \\ 1 &= z_A(t) + z_B(t) + z_C(t) \end{align*}\end{split}\]

where $p_1=4$, $p_2=2$, $p_3=40$, and $p_4=20$ are parameters with the appropriate units. $z_A(t)$, $z_B(t)$, and $z_C(t)$ are time varying concentrations of species $A$, $B$, and $C$ respectively.

Index analysis#

Determine the index of the above differential algebraic equation (DAE) system above.

Tip: do this on paper. On Gradescope, there will be a separate assignment for you to turn in your handwritten work.

Model reformulation#

Apply the index reduction algorithm from class as needed. Ultimately identify two versions of the model: one that is index 1 and another that is index 0. Find a consistent initial condition $z_C(0)$.

Implement index 1 model in Pyomo#

We will be building a library of functions.

Create model and set initial conditions#

def create_model():
    ''' Create index 1 model and set initial conditions
    
    Return:
        m: Pyomo model
    '''
    
    # Tip: Set time to go from 0 to 1 when creating the model.
    
    m = pyo.ConcreteModel()

    m.t = ContinuousSet(bounds=(0.0, 1))
    
    # Add your solution here
    
    return m

Simulate, discretize, and initialize collocation model#

def simulate_discretize_model(m,NFE,initialize):
    ''' Simulation, discretize, and initialize the Pyomo model
    
    Arguments:
        m: Pyomo model
        NFE: number of finite elements to consider (integer)
        initialize: if True, initialize the discretized model with the 
             integrator solution (boolean)
    
    Returns:
        sim: Simulator object from Pyomo.DAE
        tsim: Timesteps returned from simulator
        profiles: Results returned from simulator
    
    Overall Steps:
    1. Create Pyomo.DAE simulator and integrate with casadi/idas
    2. Transform model using 'dae.collocation' strategy. Use 3 collocation points
        per finite elemebt
    3. If initialize is true, call 'sim.initialize_model()'. This will use the
        Simulator solution to initialize the discretized Pyomo model. Really cool!
    '''
    
    # Add your solution here
    
    return sim, tsim, profiles

Plot results#

def plot_result(m, sim, tsim, profiles, include_model_values):
    """ Plot the results from the simulator (and optionally Pyomo model)
    
    Arguments:
        m: Pyomo model
        sim: Pyomo.DAE simulator
        tsim: timesteps from simulator
        profiles: results from simulation
        include_model_values: if True, also plot the values from the Pyomo model m
    
    Returns:
        nothing
        
    Actions/Steps/Tips:
    1. Plot the results stored in tsim and profiles as solid lines. Recycle code from class.
    2. If 'include_model_values' is true, plot za, zb, and zc values stored in Pyomo model 'm'.
        Use a solid symbol.
    3. Add a legend and axes labels
    """
    
    # Add your solution here
    
    # Tip: Do not forget to include `plt.show()` (assuming you imported matplotlib.pyplot as plt)

def solve_model(m):
    """ Solve discretized model with Ipopt
    
    Arguments:
        m: Pyomo model
        
    Returns:
        nothing
    """
    
    # Specify initial conditions
    def _init(m):
        yield m.za[0] == 1
        yield m.zb[0] == 0
    model.initcon = pyo.ConstraintList(rule=_init)

    # Solve collocation formulation (no objective, we are just simulating)
    solver = pyo.SolverFactory('ipopt')
    solver.solve(model,tee=True)

Simulate and solve Pyomo model with initialization#

# Create Pyomo model
model = create_model()

# Initialize discretized model with simulation result?
init = True

# Number of finite elements
NFE = 6

# Simulate model
sim, tsim, profiles = simulate_discretize_model(model,NFE,init)

# Plot simulation results
plot_result(model, sim, tsim, profiles,True)

# Solve collocation formulation with Ipopt
solve_model(model)

# Plot results
plot_result(model, sim, tsim, profiles, True)

Simulate and solve Pyomo model without initialization#

Repeat the analysis from above, but do NOT initialize the discretized model with the simulation result. (Tip: you only need to change two small things.)

# Add your solution here

Discussion: Does initialization matter?#

Write 1 to 3 sentences for each of the following questions:

What happens if we disable initialization? Does the number of interactions Ipopt need change?

Why is this specific problem robust to poor initialization? What about this specific DAE system makes it easy to solve?

Degree of Freedom Analysis#

Please answer on paper and turn in via Gradescope.

Rerun the code above for two different numbers of finite elements. Record the total number of variables and equality constraints.

How many degrees of freedom are in the simulation problem? (1 sentence is fine.)

Choose $N=3$ or a simular number of finite elements. Using the collocation equations from class, show that the discretized Pyomo model has the correct number of variables and algebriac equations. This will be a little tedious, but good to practice once on a simple model. You might need to do this when debugging a DAE model for research.

Note: You might get a strange answer. That is okay.