Steady state: The next big thing?
How to use Modelica for steady-state applications.
Modelica is a great technology for dynamic system simulation. Its tools are strong in simulating dynamics because they efficiently solve hybrid differential algebraic equation (DAE) systems. But people now want more than this.
For a number of years, I’ve heard more and more people talking about steady state, and how important it is. At the same time, people are having difficulties converging such problems. So, why are people into steady state, and why is it so difficult to compute?
The steady-state attraction
If you are a controls specialist, then you may wonder whether there is anything useful to do with a steady-state solution beyond initializing a dynamic model. If you are a system architect, however, you may be interested in steady state for design space exploration. The reason is that transients blur the picture for systems that are sized based on static performance.
More and more system architects are looking at Modelica-based applications. The reasons are simple, and the benefits are the same as for dynamic system simulation.
Modelica successfully promotes reuse, it enables highly efficient code generation, and it allows the user to work on a high level of abstraction (think equations instead of algorithms). Additionally, it is an open standard. This is pretty attractive in comparison to some legacy steady-state codes, independent of whether they were written 40 years ago in Fortran or 10 years ago in Matlab.
Robustness problem
But there’s a catch: robustness. Steady-state solutions are difficult to compute for a simple reason that is intrinsic in the problem statement.
In dynamic simulation, the state variables are known from the start and evolve along with the governing equations, whether it’s an ordinary differential equation or continuous DAE. Simulation may be inefficient, and there may be numerous challenges in detail, but in general shortening step size is a straight-forward means to achieve robustness.
All classic algorithms for steady-state simulation require a guess of the solution, and the algorithm iteratively improves the vector of unknowns. For an arbitrarily nonlinear problem, the shape of the residuals may change unpredictably a few units away from the current iterate, and the solution may be right behind the next bump or miles away.
For me, the key question to support a trend towards more steady-state applications in Modelica is how to deliver a powerful steady-state capability.
Modelica solutions
In future blog posts, I’ll comment on what tool vendors provide, and suggest future directions. For the moment, however, let’s look at what we can do with today’s Modelica tools to solve steady-state problems.
- Formulate the problem in algebraic equations
DAEs cover algebraic equations. You can thus simply skip the Modelica derivative operator der() altogether and pose your steady-state problem in its natural form. The Modelica compiler then has to apply its non-linear algebraic equation solver. All Modelica compilers I am aware of use one of the classic algorithms mentioned above. As they all require a guess of the solution (the “start iterate”), I suggest you implement a smart way to set these in your model. This works pretty well for most small problems, and converges quickly. The drawback is that in case of failure there is little to do but fiddle with the start iterate.
- Relaxation transient
A few paragraphs before I wrote about how robust dynamic simulation is. Why not exploit this for steady state and just simulate your system until the transients are damped out of the system and reach steady state? To do this, just supply an approximation of the solution values for the state variables and simulate. In theory, this actually works for both static and oscillating steady state, but takes forever. This is the show stopper for the approach. If you want to do something useful with your model, you will not have the time to wait until the relaxation transient converges. Additionally, if your approximation of the solution is not that good, you may prescribe fairly substantial transients, and your integrator starts to break down, making the whole procedure even more inefficient.
- Pseudo dynamics
The first two approaches are often combined. To do this for your model, look for the key coupling variables among the algebraic equations. For most problems, a few key variables link these equations together. If you are unsure about the physics of your problem, start looking at the set of variables selected by the tearing algorithm of your Modelica compiler. Variables of that set can be a starting point. Then, write a simple first order differential equation into your component models for each key coupling variable y, e.g., timeConstant*der(y) = y_reference – y. The correct reference value y_reference can be computed from an implicit or explicit equation system. Set a meaningful initial value for y, and see how you reduced the size of the algebraic equation system substantially.
Next time I’ll blog about what types of algorithms domain-specific simulation software packages apply to deliver a robust steady-state simulation capability.