Introduction
Reduced order models are computationally light models derived from more complex physical models. The goal of this reduction is to decrease the time needed to simulate the models. This is typically done by solving a given problem multiple times under a set of parameters, and then deducing from that a simpler model.
ROMs may have different usages:
- They can be useful if real-time computation is eventually needed. For example, when building a digital-twin, people need to simulate in real time the functioning of the device. This cannot be done through normal simulation techniques, so building a ROM is relevant in this case.
- When building a very large, system level simulation, it can be useful to simulate each piece and build ROMs to then assemble them together. This will dramatically lower the simulation time and memory usage of the entire system’s model.
- Optimizing the design of a piece typically requires performing a lot of simulations to fine tune all the parameters. When doing this, it can be relevant to use the results from the first simulations to accelerate future ones.
As you can already imagine, since this technique involves performing a number of training simulations, distributed computing on the cloud is very relevant in that case.
In this article, we will see how the main reduction technique works. Then we will show a concrete application of this technique and an implementation based on Python and the libraries PyMOR and FeniCS, with calculations done on Qarnot’s platform.
Mathematical background
In this section, we will present the main reduction technique. This will require explaining how physical equations are typically solved through numerical techniques. Don’t worry if you are not familiar with math, we will not dive into details here, but only present the ideas used.
Partial derivative equation
First, let’s define the kind of problem we are trying to solve. Physical phenomena such as the spread of heat through a metallic part, the flow of air around a plane or the propagation of electromagnetic waves are all mathematically described by equations. These state the relation between the variations in time and space of physical quantities. For example, the equation of heat tells us how temperature changes in time according to its spacial distribution at one point of time. In three dimensions, it can be written:
Discretization
There are two things to keep in mind about these equations:
- They are equations on functions. The unknown is a function that associates a value to every point in a domain and every time in a range. For example, if we are to look at the spread of heat, we would be looking for a temperature function T(t, x, y, z). These problems thus involve an infinite quantity of data. Conversely, we are not looking at models that only describe macro-quantities like the speed of the whole plane. This infinite nature makes the equation quite difficult.
- These equations link the variation of temperature in time with its variation in space. That’s what makes them partial derivative equations (PDE). To understand why these links between partial derivatives appear, let’s take the example of temperature. The variation of temperature in time (e.g. its time derivative, ∂T/∂t) depends on the flow of heat. Yet, heat flows from high temperature points to low temperature points. Therefore, the expression of heat flow involves the variation of temperature in space (∂T/∂x, ∂T/∂y, ∂T/∂z). The fact that the equation links different partial derivatives together makes it all the more difficult.
It is generally impossible to mathematically find the solution of such a partial derivative equation. In fact, some common PDEs are so hard that mathematicians struggle to demonstrate anything about them. This is the case of the Navier-Stokes equation (describing the flow of a fluid) where it has not been shown yet that a smooth solution always exists (a prize of 1 million is actually on this problem).
Therefore, engineers and scientists look for approximate solutions through numerical means. In these problems, some parameters may come into play and can change. They can be for example physical property of the material that can change if we decide to replace one material by another or if we study the use of different possibilities. Parameters can also be the ambient temperature, or the plane speed, etc. For these reasons, we may have the same geometry and the same model, but we would need to compute the solutions every time we want to study what happens for a new set of parameters.
To solve equations numerically, the first step is always to derive a finite model from the infinite equation given. Indeed, as we have seen, the models we are trying to solve contains an infinite number of unknowns since we must find the values at each point and time. This can’t be done by a computer, since it only has finite memory and computational power. Therefore, we need to make this problem finite. Several techniques can be used to perform this discretization process, but they often share the same two steps: meshing and equation discretization. Both these steps are explained underneath.
Meshing
The idea of this step is to say that since we cannot calculate the values at every point and time, we will only compute it at some points and time values. Then, when we need to know the value at a given point, we will infer it using the neighbor points where the values have been calculated. For example, it can be calculated through linear interpolation (e.g. taking the balanced mean of neighboring points).
A basic example is as follows. Consider a thin bar of heated metal of length 10. Then, the piece has only one dimension, and we could divide the length in a number of segments (let’s say 5). If we know the temperature at every one of the extremity points, we can create a linear curve that will be close to the real solution.
The central idea in this process is that now, the solutions we are looking for are fully described by a finite number of information. Here, it is the values at 0, 2, 4, 6, 8 and 10. Thus, we can say that the space we are searching a solution in now has a finite dimension (6 here), or that the solution has a finite number of degrees of freedom (DOFs).
In higher dimensions, discretizing space often involves building a mesh, which is a partition of the object into small triangles or tetrahedrons. On each of these primitives, the final function will be linear. Here is the meshing of a 3D piece.
On the other hand, time is very easily discretized by cutting the total time windows into small time steps.
To understand reduction order techniques, it is important to note that it is analogous to say that the function is an interpolation of values at some points and that it is a linear combination of some simple functions. If we take back our example, we can see that the total solution is the sum of the six function represented in dots:
Every one of this basics functions are the scaled version of a more “primitive” functions whose values are 1 at a point and 0 at each other point (here the value is 10, so we can see better).
In the mathematical languages, these 6 basic functions are called a basis of our solution space. This means that we are looking for solutions that can be expressed as linear combinations (e.g. a sum of these functions each multiplied by a scaling factor) of these functions.
Equation derivation
Now that our problem has been restrained to a limited number of degrees of freedom, we must adapt the equation to that sort of solutions. There is a principle in math that states that for a system to have a unique solution, there has to be as many equations as degrees of freedom.
There are different techniques that can be used to do so. The simplest one would be to rewrite the equation at each point, approximating space derivative by (T(x+Δx) — T(x))/Δx, e.g. by the difference with neighbor points. This is known as the finite difference method and though it is the simplest one, it is only rarely used. The most used one are the finite element method (FEM) and the finite volume method (mainly for fluid dynamics).
Below is a quick explanation of the finite element method. If you are not interested in the details, or you already know about it, you may jump to the next section that explains the core idea of model order reduction.
The finite element method
To understand the idea behind the finite element method, we can go back to our heat equation for a metal bar and the more natural finite difference technique. In this technique, we approximate the derivative in one point by the difference with the neighboring points : ∂T/∂x = (T(x+Δx) — T(x))/Δx where Δx is the distance with the next point. The problem is that, as we see, we divide our difference by Δx to have a slop. If one wants a fine mesh for a precise calculation, then Δx will be very small. Because we divide by that very small number, ∂T/∂x can get very big and very sensible to numerical error. In the heat equation, it is even worse because we consider the second derivative, which is in Δx².
Therefore, the finite difference method is quite instable. One of its other disadvantages is that there are not many theorems on that technique. For example, we have little results about the error we make when using this technique.
To tackle the instability issue of the finite difference, the idea of the finite element method is to use the opposite operation of derivation, integration. Since integration operates like a mean, it seems logical that it makes for more stable models. However, integrating over the whole domains makes for a loss of information. Indeed, it is not because two function have the same means that they are equal. To solve this and preserve the richness of our function equation, we multiply both sides by a function, called the test function (written v in most cases). It can be seen as a windowing function that will only preserve a part of the unknown function. This time, if two equation have the same mean over each window, then we can expect them to be equal. In fact, under the right hypothesis, we have an equivalence between the two formulations. For example, considering the heat equation in steady state, the two following formulation are equivalent (where V is the appropriate space of test function):
The equation on the right-hand side is called the weak form equation and is the one that is used in the finite element method. After modifying the expression, the key idea of FEM is not only to search for a solution in a space Vh of finite dimension n, but also to test it with function v in that same space. Since Vh is of finite dimension and the equation linear in v, verifying the equation for all v in Vh is equivalent to verifying it for all element of its basis (v1, …, vn). That leads naturally to a set of n equations. By considering the n coefficients λ1, …, λn so that T = Σ λivi, we obtain a system of n equation with n unknowns. This system can be systematically solved.
As a conclusion, the finite element method consists in rewriting the equation in another form, the weak form, which involves a test function. The discretization is done by using test functions that have the same degrees of liberty as the solution space. For mathematical reasons, it happens that, because we have used a weak formulation, we have a lot more mathematical results and a much better understanding of that technique. For example, Céa’s lemma gives an upper limit for the error we make using this discretization technique. More precisely, it gives an upper bound for the distance between the real solution and the best approximation of that solution in our solution space.
Model order reduction
The core idea of model order reduction is to look for a solution in a specific basis that contains only a few elements but that we know will suffice to describe the solutions. The type of basis exemplified in the meshing section is very relevant when we don’t know what the solution will look like. Indeed, it offers a good flexibility and we are sure that we will be able to describe well any function on the domain. Yet it is also quite naive:
- First, we can expect the degree of liberty of the solution to be somewhat comparable to the number of parameters. A mesh can contain millions of points, while we often consider only a handful of parameters.
- Second, the elements of that basis are simple and convenient for calculation but not physical at all. Indeed, in these elements, values at a point are totally uncorrelated to the neighboring points while in reality there is strong link between values at points that are near. To give an example, the same thing happens with numerical images. The most basic way to describe one is to tell the color of each pixel. It is a convenient and simple description, but it doesn’t use the fact that the color of one pixel is typically close to the ones around. From that assessment, people had the idea to express an image in some other basis (using Fourier decomposition) to compress the size of the image. This is the core of image compressing and file format like .jpeg
Thus, the idea is to simulate some solutions using a full-order model and then to deduce from it a smaller basis that will be able to express solutions close to the full order ones. Then, we will only need to derive equations in this new solution space. Since this space has a limited number of degrees of freedom, the derived model will only consist in that same low number of equations. Therefore, it will be almost immediate to solve. Moreover, since the solution space was built to be able to describe solution that are close to the real one, we know that the solution derived from our reduced order model will be close to the full order one.
There is not a lot more to know about the basics of model order reduction. We may now explore some ways to build the reduction basis. The most basic way to build it is simply to do some simulations and use directly the solutions to make out the basis (we would only orthonormalize the basis). However, it is possible to do better. Below is a quick description of some algorithms that can be used.
Some algorithms to build the reduction basis
Proper Orthogonal Decomposition
This technique is close to a technique called PCA (principal component analysis), which is a very common technique in data science. The idea is to extract a basis of lower dimension that describes as well as possible some vectors. So, when using POD, the pipeline is first to perform a number of simulation (let’s say 100) and then to extract a smaller basis (say 20) that describes as well as possible the results obtained.
Strong greedy algorithm
In this algorithm, the idea is to build the basis steps by steps. Starting from a set of results, we will choose one that will be our first vector of the reduction basis. The reduced order model related to the basis is created and one can compute the solutions for all the other values from the set of parameters. Then, the solution with the biggest difference between full order and reduced order is added to the basis. In doing so, we lower the error for the solution that was the least well described. This operation is then repeated over and over again till we reach the desired number of vectors.
Weak greedy algorithm
The two aforementioned algorithms required computing the solutions for every parameter of a training set. In the weak greedy algorithm, we will follow the approach of the strong greedy algorithm, yet, we won’t need to compute all the exact solutions because the error is not computed but solely estimated from the ROM.
To learn more about these algorithms, I recommend reading this tutorial, which goes deeper into how these works and their advantages and drawbacks.
Going further
If you are interested in implementing Reduced Order Model computations in the cloud, have a look at our blog or at Qarnot website.
