Case study 2 – Solving algebraic equations (multi unit)

Key concepts

  • Solving nonlinear algebraic equations
  • Systems of equations
  • Degrees of freedom
  • Slicing lists
  • Arrays
  • fsolve


When analysing plant data, rarely are we provided with all rate and composition data. Instrumentation, sampling and laboratory analysis can be expensive and we need to make best use of available information to infer data for which we don’t have measurements. In this post we will look at an example where we exactly the minimum amount of the information to complete the mass balance across a system of 3 separation columns. In a later post, we will look at managing situation when we have either too many or too few measurements to complete the mass balance.

Problem description

The problem for this post is the 2nd example taken from the collection originally distributed at the ASEE Chemical Engineering Summer School held in Snowbird, Utah on August 13, 1997. This problem involves an array of 3 distillation columns separating xylene, styrene, toluene and benzene. The arrangement of the columns and the names for all the streams are shown in the following diagram.

Distillation column array

We are given the molar flowrate for stream F (=70 mol/min) and molar fraction for all 4 components in streams F, D1, B1, D2 and B2.


We are then asked to provide the following information:

  1. Molar flowrates for product streams
  2. Compositions and molar flowrates of the intermediate streams

It is possible to solve this problem separately in these 2 steps. However, in this blog we will create a single system description which will provide molar flowrates and compositions for all streams at the same time.

Degrees of Freedom

Understanding the degrees of freedom is critical for a problem like this. It may not be immediately obvious wether or not we have enough, too little or too much information. In fact in this case the problem is over specified by 5 results, but I am getting ahead of myself.

First of all, how many variables do we have? For each stream there are 5 variables = 4 molar fractions + 1 molar flowrate. Also, there are 7 streams = 1 feed + 2 intermediate + 4 product streams. So in total there are 35 variables.

How many equations do we have? Before describing the equations, we can use the following logic. For each distillation column, there are 4 mass balances, corresponding to each of xylene, styrene, toluene and benzene (note – this purely a physical separation and there are no chemical reactions). So over the 3 columns, there are 12 mass balance equations. In addition we can assume that there are no other chemicals involved so that the molar fractions for each stream must sum to 1. This creates an additional 7 equations.

Finally, as described above we have been given 21 pieces of information = 20 molar fractions and 1 molar flowrate.

So our degrees of freedom = 35 – (12 + 7) – 21 = -5. Or in other words it has been over specified by 5 variables. This is because we are given all molar fractions for each of the 5 streams, when we could in theory calculate the fourth fraction given only 3. So the problem is only trivially over specified. In otherwords, we can only use the additional 5 data points to check that the supplied molar fractions sum to 1. As they do, we really have a exactly specified problem.

Describing the mass balance

In this blog, we will deviation from the problem description in the original paper. Instead, we will take advantage of python’s ability to use vectors, to create a simpler, more generic and frankly more robust description.

F * F_{mf} = D * D_{mf} + B * D_{mf}

\sum F_{mf} = 1.0

\sum D_{mf} = 1.0

\sum B_{mf} = 1.0

  • F, D, B = molar flowrates, mol/min, for streams F, D and B
  • F_mf, D_mf, B_mf = vector of mole fractions for streams F, D and B

The first equation represents 4 different equations, corresponding to the 4 chemical components. Also, this set of equations is repeated for each of the 3 distillation columns.

What type of problem

As with the case study 1 (previous blog post) we have a number of nonlinear algebraic equations to be solved. Algebraic because there are no differential terms. In this case the the equations are nonlinear because they involve multiplication of different variables for which we are seeking a solution (e.g D * D_{mf}).

Knowing this, our general approach for solving will be:

  • Define known variables
  • Create a function containing all equations in difference form (i.e. LSH – RHS)
  • Solve using a routine for finding the roots of equations (in this case fsolve)

The solution

To begin with here is the full solution, produced using python through a Jupyter notebook.

Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.


As with all of our python solutions, we begin by specifying the individual packages we will need.

import numpy as np
from scipy.optimize import fsolve

This time we will use the following familiar packages (also used in the post for Case 1).

  • numpy – a very powerful library of mathmatical routines
  • scipy – an extensive library of mathematical analysis routines

We then continue by specifying the known molar flowrate for stream F and mole fractions for the feed stream and 4 product streams.

F = 70.0
F_mf  = np.array([0.15, 0.25, 0.40, 0.20])
D1_mf = np.array([0.07, 0.04, 0.54, 0.35])
B1_mf = np.array([0.18, 0.24, 0.42, 0.16])
D2_mf = np.array([0.15, 0.10, 0.54, 0.21])
B2_mf = np.array([0.24, 0.65, 0.10, 0.01])

As mentioned above we are going to use vectors to represent the set of mole fractions for each stream. This is done first of all by creating a python list – i.e. a list of items separated by comma and contained in square brackets. This list is then converted to an array using the numpy function called array. However, rather than doing this as 2 steps, this is all done in a single line of code. The numbers in the code should correspond to the table in the original problem specification.

Finally we need to specify initial guesses for all the remaining variables.

D,B,D1,B1,D2,B2 = [10.0]*6
D_mf = np.array([0.1]*4)
B_mf = np.array([0.1]*4)

The code takes advantage of some of the syntax shortcuts available in python. In this case, all of the unknown molar flowrates are initialised to 10.0 mol/min in a single line of code (note: the code [10.0]*6 creates a list containing 6 values of 10.0 = [10.0, 10.0, 10.0, 10.0, 10.0, 10.0]).

Similarly the mole fractions for the intermediate streams are all initialised to 0.1 and put into an array.

The mass balance

The heart of the solution is contained in the next bit of code, which describes the model equations:

    error_1 = 1 - sum(D_mf)
    error_2 = 1 - sum(B_mf)
    error_mf_1 = D *D_mf   + B *B_mf   - F*F_mf
    error_mf_2 = D1*D1_mf  + B1*B1_mf  - D*D_mf
    error_mf_3 = D2*D2_mf  + B2*B2_mf  - B*B_mf

These equations are the mass balances for the 3 distillation columns and the specification that the mole fractions must sum to 1.0 for the unspecified streams. Also the equations are written as a difference between the left hand side and right hand side that they can be solved through numerical solution.

These equations need to contained in function which will be called by fsolve. We have 14 variables which need to be solved for and they all need to be wrapped into a single list to be passed by the fsolve function. To do this we create a function called mass_balance which accepts an argument x. It is x which is a list containing the 14 variables. However, we need to split x apart so that it fits into out mass balance equations.

def mass_balance(x):
    D, B, D1, B1, D2, B2 = x[0:6]
    D_mf = np.array(x[6:10])
    B_mf = np.array(x[10:15])
    error_1 = 1 - sum(D_mf)
    error_2 = 1 - sum(B_mf)
    error_mf_1 = D *D_mf   + B *B_mf   - F*F_mf
    error_mf_2 = D1*D1_mf  + B1*B1_mf  - D*D_mf
    error_mf_3 = D2*D2_mf  + B2*B2_mf  - B*B_mf
    return np.append([error_1, error_2], [error_mf_1, error_mf_2, error_mf_3])

You can see that this is done by using the notation x[start_position, final_position]. However, this highlights 2 conventions that you need to get used to with python:

  • The first position in python is always referenced by 0 (not 1)
  • In sequences the final position is not included

So x[0:6] actually slices out the first 6 elements of the list x. That is x[0], x[1], x[2], x[3], x[4] and x[5]. Equivalently, x[6:10] slices out 4 elements starting from the 7th position. That is x[6], x[7], x[8] and x[9].

Depending on you background this may at first be confusing. Take some time to remember this as it is used throughout python and is simple once you get used to it!

Finally, the function returns the 14 unknown variables grouped together into a single list. The order of this list does not have to relate to the order of variables contained in x, but there does need to be the same number (i.e. the degrees of freedom).

Solving the equations

Now that the problem has been properly specified, the code to solve the problem is relatively simple as it uses the power of the fsolve routine.

solution = fsolve(mass_balance, np.append([D, B, D1, B1, D2, B2],[D_mf, B_mf]))

D, B, D1, B1, D2, B2 = solution[0:6]
D_mf = solution[6:10]
B_mf = solution[10:15]

This code calls the routine fsolve and tells it to use the model contained in our function, mass_balance and uses the initial values we set for the molar flowrates and mole fractions. It returns the numerical result as a list to the variable solution. From here we slice out the values for the molar flowrates and mole fractions.

Finally we are able to display the results. It is always worth sense checking the results. In this case a quick check of the following all works out:

  • Flows in and out of each column are equal
  • Mole fractions sum to 1.0 on each stream
  • Random selection of chemicals in and out of a column (e.g. xylene in and out of column 1) are equal

Sense checking and interpretation

Let’s do a quick sense check on the results to ensure that there is nothing obviously wrong with our solution.

  • Results are all positive (check)
  • Molar flowrates sum across each column (check)
  • Mole fractions sum to 1.0 for each stream (check)
  • Random selection of individual mass balances are correct (check)
  • Flowrates are matching the reporting location of different chemicals – e.g. F is high in toluene, D is also high in toluene and the majority of the florwate goes to D as expected (check)

Case Study 1 – Solving algebraic equations

Key Concepts

  • Solving nonlinear algebraic equations
  • Degrees of freedom
  • Python packages
  • Functions
  • fsolve
  • Iterating lists


Solving a nonlinear equation is relatively common task for process and chemical engineers not only in technical roles, but also for supervisory, operating and strategic roles. I’ve seen people use the “guessing” method to hone in an answer and others pull out text books to remember how to find a algebraic solution to the equation. Finding numerical solutions to nonlinear equations should become second nature and we will look at how to do it here in a way that clearly lays out the solution so that you and others can review it later on.

Problem description

The problem used in this post has been taken from a collection (1st problem) originally distributed at the ASEE Chemical Engineering Summer School held in Snowbird, Utah on August 13, 1997. The problem requires calculation of different gas properties based on the van der Waals equation of state:

The problem requires us to solve for molar volume, V and you can see from this equation that it is not trivial to rearrange this equation and solve for V. Instead we will use a numerical equation solver to find a solution.

Before progressing to the code, we need to ensure that we have a properly defined problem. The van der Waals equation contains 6 variables:

  • P = pressure (atm)
  • V = molar volume (liters/g.mol)
  • T = temperature (K)
  • R = gas constant (atm litre/g.mol K)
  • a, b = parameters to the equation

We are also given the following equations:

a = \frac{27}{64} (\frac{R^2 T_c^2}{P_c})

b = \frac{R T_c}{8 P_c}

These equations introduce 2 additional variables:

  • Pc = critical pressure (atm)
  • Tc = critical temperature (K)

Degrees of freedom

In total we 8 variables and 3 equations, leaving us with 5 degrees of freedom. In other words we need define 5 further pieces of information to solve for V. In this case we are given the gas constant R (or we could have looked it up), we are told the gas is ammonia and actually given Pc and TC (or we could have looked them up) and for each part of the problem we are told the conditions T and P (or in some cases we are given the reduced pressure Pr).

So now we have a fully defined problem and are able to calculate V. However, we will also be required to calculate the compressibility factor, Z using the following equation:

Z = \frac{P V}{R T}

This equation only introduces the new variable Z. This is 1 new equation with 1 new variable so there are zero degrees of the freedom and the problem remains fully defined.

What type of problem

Before continuing we need to think about what type of analytical problem do we have here as this will help us understand the core concepts required to find a solution, plus has implications for how to setup with solution approach in Python.

In this case, it is fairly simple as we have a number of nonlinear algebraic equations to be solved. Algebraic because there are no differential terms and nonlinear because the term V (i.e. the variable we need to be solving for) appears with a squared term and we can not rearrange the equations with V having only a power of 1.

So our general solution approach will involve:

  • Define known variables
  • Create a function containing all equations in difference form (i.e. LSH – RHS)
  • Solve using a routine for finding the roots of equations (in this case fsolve)

The solution

I will show the entire solution to the problem before explaining some of the fundamental solution concepts. This output has been produced using python through a Jupyter notebook.

Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

Python packages

The first concept introduced in the solution is the use of packages. Python is program language which ships with a well developed set of commands and functions. These provide a lot of functionality. But the real power of Python comes from the ability to pull in packages to perform an ever expanding range of functions. Fortunately, the python packages for data analytics are extensive and very mature.

The first block of code imports from 3 packages which get used extensively in data analytics:

  • numpy – a very powerful library of mathmatical routines
  • scipy – an extensive library of mathematical analysis routines
  • matplotlib – a flexible charting and visualisation library
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import numpy as np

I will right a separate post about importing libraries, but the key point to understand in this example is that:

  • A procedure called fsolve is now available to the code
  • The pyplot part of the larger matplotlib library is now available to the code and can be referenced through the prefix plt
  • The numpy library is now available to the code and can be referenced through the prefix np

Defining functions

The next part of the code simply defines the supplied constants. Then we get to the following step:

def vdw_calc(V):

    a = 27.0 / 64.0 *(R**2 * Tc**2 / Pc)

    b = R*Tc / (8.0 * Pc)

    converg_error = (P + a/V**2)*(V-b) - R*T

    return converg_error

This code creates a function called vdw_calc which accepts a value for the molar volume V. It then uses the equations for a and b, utilising the values for all the parameters which will be defined before this function is called for the first time.

The actual van der Waals equation of state is rearranged so that the right hand side (RHS) is subtracted from the left hand side (LHS). The difference is called converg_error and our goal will be to find a solution which makes converg_error = 0. When this occurs, then we have a solution to the van der Waals equation. If the value for V supplied to the routine is not an exact solution, the converg_error represents how far off the van der Waals equation is from being solved. This value is returned and will be used later in the code.

Similarly a procedure is created for calculating the compressibility factor. This procedure is called comp_factor and does not receiving any values, instead using the pre-defined values for P, V, R and T to calculate the compressibility factor, Z. The value for Z is returned.

def comp_factor():

    Z = P*V / (R*T)

    return Z

Solving part a) – Solving the nonlinear equation

Part a) of the original problem was simply to calculate the molar volume, V and compressibility factor, Z at a specified T and P. The bit of code that does this is:

T = 450.0

P = 56.0

V = fsolve(vdw_calc, 1.0)[0]

Z = comp_factor()

The heart of this solution is using fslove (from scipy library) in combination with vdw_calc (our routine which takes a value for V and returns the error between the LHS and RHS of the van der Waals equation). fsolve accepts 2 bits of information. Firstly it needs to know the equation it needs to solve for (in this case vdw_calc) and a starting approximation for the solution (in this case we are guessing that the right answer for V is 1.0). fsolve then iterates different values for V (starting with 1.0) until vdw_calc returns 0.0.

Once fsolve has a solution it returns the correct value for V and assigns it to our parameter V. fsolve actually returns the answer as a list, just in this case it is a list with only a single number in it. We want V to contain just the number V and not the list. We extract the number by using the syntax [0], which literally means give me the first number in the list (python numbering always starts with 0 not 1).

Now that we have V and using the same values for T, P and R we calculate Z. The last 2 lines of code simply display the results with a little bit of formatting. I will not explain how this works for the moment,

Solving part b) – Solving the equation multiple times

Part b) of the original problem requires essentially the same calculation to be repeated for different values of pressure, P. In this case the pressures are given as reduced pressures, Pr and will need to converted to pressure according to the following equation:

P = P_r * P_c

Python provides a simple construct for storing information called a list. In this case we can store the list of reduced pressures we need to solve for in a list called Pr_list. Once this is defined, we can iterate through each value using a for loop:

Pr_list = [1,2,4,10,20]

for Pr in Pr_list:

    P = Pr * Pc

    V = fsolve(vdw_calc, 1.0)[0]

    Z = comp_factor()

In this case, the variable Pr takes on each value, one by one in the list and is then used to calculate the pressure, P, which we need for the van der Waals equation. There is no need to understand how many items are in Pr_list or how many times the for _ in loop has repeated. Python handles all of this seamlessly.

The rest of the code is essentially the same as part a), except it is repeated for each iteration.

Sense checking and interpretation

When ever solving analytical problems, always, always sense check your answer. There are many techniques you can employ depending on the situation. For example, in this case we can check the following:

  • Results are all positive (check)
  • Molar volumes are right order of magnitude (check)
  • Compressibility factors are right order of magnitude (well actually I don’t know what to expect)
  • Nonlinear form of the compressibility factor with pressure is expected (the term \frac{1}{V^2} is consistent with high volumes at very low pressure – check)
  • Manual check of some results using a calculator (check)

There is not a great deal to interpret here with out further context. However, I will revisit this problem in the future with opportunities to scale the approach for a more industrial environment.

Predictive Maintenance Obstacles

Predictive maintenance is an important technique which can and should be a fundamental capability for engineers working industry. Unfortunately few engineers are aware of the techniques and even fewer understand the tools & techniques.

The mathematics of remaining useful like (RUL) estimation is reachable for a wide range of engineers and I encourage everyone to increase their awareness.

Mathworks Predictive Maintenance Fig2

Wether implementing your own predictive maintenance solutions or working with others, you need to be aware of the typical obstacles to achieving a useful outcome. Four  major obstacles are raised in the post, Top Obstacles to Overcome when Implementing Predictive Maintenance.

These are:

  1. Being Unaware of How to Do Predictive Maintenance
  2. Lacking Data to Create Proper Predictive Maintenance Systems
  3. Lacking Failure Data to Achieve Accuracy
  4. Understanding Failures but Not Being Able to Predict Them


Full Potential and Latent Capability

Understanding the full potential of your system will lead to opportunities in latent capability. Opportunities that are generally low cost and low risk, as they take advantage of assets you have already installed and have a history with operating.

However, understanding how to release this capability and being confident in what a system is capable of requires a range of different skills. For example how can you turn the noise of the following data into information that can be used change your business?

full potential time series

One technique that can be used with limited data is a Weibull plot of operating data……

full potential weibull

I will be sharing insights on a range of techniques (including Weibull) for understanding full potential during this FREE IChemE sponsored Webinar .

Attendees will get to appreciate why understanding the Full Potential of a process as a whole is so important and techniques for analysing operations and integrated value chains. This is a critical capability for businesses wanting to understand the value of untapped potential and how it links to targeted improvement activity. 

Webinar – Understanding Full Potential





Understanding Full Potential (IChemE Webinar)

“Latent Capacity” aka untapped potential is ubiquitous across the mining and mineral processing industry. We know it’s there, but we can’t always see it. When we see, we don’t always know where to start. And when we start, we find that change is a difficult road. This untapped potential is not due to a lack of effort. Operations and providers to the industry alike are focused on improvement and “closing the gap”. For many, Big data and advanced analytics provide hope for tantalising returns from existing assets. While everyone agrees the stakes are high, the pathway to realising the full potential is fraught with traps for the unwary hence experience is everything.

I will be sharing some insights during this FREE IChemE sponsored Webinar aimed at engineers dealing with the challenges of process optimisation.

Attendees will get to appreciate why understanding the Full Potential of a process as a whole is so important and techniques for analysing operations and integrated value chains. This is a critical capability for businesses wanting to understand the value of untapped potential and how it links to targeted improvement activity. 

Webinar – Understanding Full Potential

Valuing Innovation


Gaining traction with innovation requires a clear and compelling case for change. A case that overcomes the opportunity cost hurdle for access to time, money, resources and leadership support. A failure to do this well means that many great ideas are left undiscovered and holders of the ideas frustrated that nobody else “gets it”.

In fact, in this dynamic and uncertain environment a good case for change it not sufficient – a great case must be presented.

A great case for innovation requires a narrative that engages the full breadth of stakeholders in a picture of what the future could be. The narrative helps win support for looking at the idea in more detail while presenting the light on the hill for your supporters.

A great case for innovation will obviously contain a robust financial analysis of the expected benefit. However, this analysis needs to be comprehensive and include likely associated costs plus potential downsides in addition to the good news story. A partial analysis will lack credibility and fail to uncover the true impact to the business.

A great case also needs to consider the range of non-financial or non-measurable benefits and costs. Often an innovative idea will look good based on the financial analysis, but it will be the linkages to a broader business strategy or simply the emotion of being seen as a leader in the industry that will make for a great case.

A great case must also consider the various risks that come from the idea itself, the implementation and the potential scenarios from changes in the environment.

There are many techniques available to innovators to deal with the various elements of a great case. The right techniques depend on the type of innovation you are considering, its impact to your business and the stage of innovation between ideation and final implementation.

This is just one topic that Geoff Carter of On Track RDI and I will explore in the Practical Innovation masterclass. The masterclass will be held in Perth and is booked for Thursday 27th October from 12 to 4pm, with networking afterwards.

Places are limited, but still available. For more information please follow this link and register here.

Practical Innovation


Australians are renowned for their smart ideas, but we often fail to back them and turn them into commercial realities.”  – National Innovation and Science Agenda Report, December 2015.

Geoff Carter of On Track RDI and I have teamed up to create a masterclass which explores a practical approach to innovation. An approach which encourages smart ideas, pragmatically priortises effort on the right ideas, and leads to a minimum viable product for prototyping. We then explore how to engage the end users in turning prototypes to a bottom line change.


This is a great opportunity for people getting into innovation programs through to people already leading innovation programs who want to improve their impact.

Some of the topics we will cover include:

  • Leadership for creativity
  • Team culture for innovation
  • Idea generation and control
  • Delivery with uncertainty
  • Rapid testing and change management
  • Comprehensive value propositions

The masterclass will be held in Perth and is booked for Thursday 27th October from 12 to 4pm, with networking afterwards. For more information please follow this link.

Register Here

I am really excited about this opportunity to present with Geoff. So join Geoff and me for an engaging and insightful workshop with a group of like minded people. Places are limited so we can have interaction and participation throughout the afternoon.