Thanks for making this material available, it has helped me a lot in my journey trough operations research and confused me at the same time, mainly because there's a couple of glitches in the code, maybe they are intentional to help a careless student, but if so, I hope that this suggestion is not out of place, it's better to leave a clue around the footer of the text in the form of a comment, the glitches:

• In google colab you have to actually install the packages for pyomo and glpk it can be easily done in the preamble by adding a comment <!pip install pyomo> <!apt-get install -y -qq glpk-utils>
• The code for the model constraints in section 6.8.3 reference an undefined set feeds it should be streams this same error is found in sections 6.8.5-7.

Wrong lines

# constraints
m.cons = pyomo.ConstraintList()
for s in S:
    m.cons.add(sum(m.x[s,p] for p in P) <= streams[s]['avail'])
for p in P:
    m.cons.add(sum(m.x[s,p]*(feeds[s]['octane'] -    products[p]['octane'])       for s in S) >= 0)
    m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
    m.cons.add(sum(m.x[s,p]*(feeds[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
    m.cons.add(sum(m.x[s,p]*(feeds[s]['benzene'] -   products[p]['benzene'])      for s in S) <= 0)

Corrected Lines

    m.cons.add(sum(m.x[s,p]*(streams[s]['octane'] -    products[p]['octane'])       for s in S) >= 0)
    m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmin']**1.25) for s in S) >= 0)
    m.cons.add(sum(m.x[s,p]*(streams[s]['RVP']**1.25 - products[p]['RVPmax']**1.25) for s in S) <= 0)
    m.cons.add(sum(m.x[s,p]*(streams[s]['benzene'] -   products[p]['benzene'])      for s in S) <= 0)

Hope this helps.

Since I'm reading this as a novice, please regard this petition, may you send me a concrete representation of the matrix form for this problem, set notation is unclear to me although its advantages for generalization and problem statement, are obvious, it makes it really difficult to formulate a constraint matrix, if my mistake is far worse than it seems could you please reference some texts that help with a concrete formulation in matrix form, thank you <3.

The solver has difficulty with the stiff ODE in the exothermic CSTR example.
The solution is not the same as the textbook solution for 300K, any idea for improvement?

I tried a few options in scipy for solve_ivp, namely "Radau" and "BDF", but didn't solve the issue.