SciPy is a general-purpose scientific computing library for Python, with an optimize module for optimisation. This was used to solve the water tank reference problem. Both sequential and global solutions are presented below.
One item which was unspecified in the description of the reference problem was the details of the demand. For the purpose of this exercise the demand will be cyclic, with periods of peak demand (with a sinusoidal profile) separated by intervals of zero demand.
Solvers
The SciPy optimisation module offers a number of potential solvers. We want to use constraints, so this limits our options to:
COBYLACOBYQASLSQPandtrust-constr.
All of these solvers also allow for bounds on the decision variables.
Jacobian & Hessian
The Jacobian and Hessian play critical roles in determining how an algorithm navigates the search space to find optimal solutions, particularly in gradient-based optimization methods. The Jacobian, which represents the first-order derivatives of the objective function, gives gradients, while the Hessian, which is a matrix of second-order derivatives of the objective function, gives curvature.
For most solvers the Jacobian will determine the direction of descent along the objective function surface, while the Hessian might be used to determine the step size and refine direction of descent.
The Jacobian is quick to compute, while the Hessian can be expensive. For the global problem a function is provided for the Jacobian, otherwise it would be calculated numerically. The impact on the results either way is minimal.
Global Optimisation
One approach is to make the tank water level one of the optimisation variables. This does make the formulation of constraints and objective easier, but it also makes the dimensionality of the problem higher. So the optimisation runs slower. It’s better to simply calculate the water level where required.
I set up the sequential solution to the reference problem using SciPy in a Jupter notebook. The required packages (and a few non-essentials) are in the requirements.txt.
In addition to the objective() function I have also implemented jacobian(), which is passed separately to minimize(). In principle it’s not necessary to supply an explicit Jacobian and the solver will infer the gradients using numerical techniques. However, this does mean that the objective function is evaluated many more times. So if possible it does make sense to provide an explicit implementation of the Jacobian.
SLSQP
The SLSQP (Sequential Least Squares Quadratic Programming) solver is suited to non-linear programming problems with both equality and inequality constraints. It’s effective for constrained optimization problems involving smooth objective and constraint functions.
The top panel shows the inflow rate. The second panel shows the outflow rate (blue) and demand (black). The third panel shows the tank level.
During the first cycle, since the tank starts out only a quarter full, the inflow pump works hard (but not at its maximum rate!) to satisfy the demand and fill the tank. Towards the end of the cycle the tank is essentially empty and outflow is coming directly from inflow. The optimal inflow rate during this first phase is just high enough to ensure that the demand is satisfied for the whole cycle.
When demand subsides after the first cycle the inflow rate drops to a lower level, which is then sustained for the remaining cycles. The tank level slowly rises with time, dropping every day during the period of demand and then recovering afterwards.
Trust Regions
The results from the Trust Regions (trust-constr) solver are given below. This solver uses a trust regions to iteratively approximate the objective function.
The solution from the Trust Regions solver is similar to that from the SLSQP solver. The outflow rate matches the demand. However, the evolution of the inflow rate is nowhere near as clean, oscillating during each demand cycle with some apparently random perturbations. From an aesthetic perspective (if one is permitted to judge optimisation results on the basis of aesthetics!) this solution is not as pretty as the previous one. Furthermore, as we’ll see in the summary below, the value of the objective function for this solution is also higher, so it’s also objectively an inferior solution.
Summary
The summary table below compares the solutions obtained via the two solvers considered above.
| Optimiser | Objective | Time (s) | Timestamp |
|---|---|---|---|
The SLSQP solver is both quicker and gives a better value for the objective function. The two runs were otherwise identical, only the solver was changed. I’m sure that I could tweak the parameters for these two solvers to achieve similar results more rapidly (most likely by tweaking the convergence criteria), but for the moment I’m happy with these solutions as they will provide a useful reference for comparison with the sequential optimisation results below.
Sequential Optimisation
I set up the sequential solution to the reference problem using SciPy in a Jupter notebook. I spent a lot of time fiddling with this and probably over-engineered the implementation significantly. Here are some modules that are used by the notebook:
setup.py— problem configuration (use by both sequential and global solutions)flow.pyplot.pypump.pyandtank.py.
SLSQP
These are the results for the sequential optimisation problem using the SLSQP solver. Comparing them to the global problem results above it’s evident that there are some major differences.
The “greedy” nature of the solution is readily apparent: at the start of the first demand cycle the inflow rate increases to satisfy the demand. However, it then reaches the maximum inflow rate bound and cannot go any higher. Since the demand exceeds the inflow rate the tank level drops to zero and for a while the outflow pump is unable to completely satisfy the demand. Essentially the sequential optimisation doesn’t have the foresight to increase the inflow pump rate from the start in order to satisfy the first demand cycle. This only happens on the first cycle, but illustrates the challenges with this approach to optimisation.
The inflow rate for each demand cycle establishes a regular pattern, where it increases rapidly at the start of each cycle, hits a plateau and then drops back to zero once the tank is filled. The inflow plateau is close to the maximum inflow rate.
Trust Regions
The consequences of “greedy” optimisation are again apparent, where the outflow falls below demand in the first cycle. In contrast to the SLSQP solution, the level of the inflow plateau (except on first cycle) is around 40% of the maximum inflow rate, meaning that each inflow cycle is of longer duration.
Summary
For the sequential optimisation problem it’s not possible to quote a single value for the objective function since this is treated separately for each of the time steps. As a result it’s difficult to objectively compare the quality of the solutions. However, it was apparent that the Trust Regions solver was substantially (an order of magnitude) slower.
Conclusion
Based on the results above I can see that the sequential optimisation approach is not ideal for this problem. Following posts will use CVXPY and Pyomo and will focus exclusively on the global optimisation approach.