{'fraction_non_alphanumeric': 0.04017292281737663, 'fraction_numerical': 0.012020244622522143, 'mean_word_length': 4.359988697372139, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shape-Based Nonlinear Model Reduction for 1D Conservation Laws

We present a novel method for model reduction of one-dimensional conservation law to the dynamics of the parameters describing the approximate shape of the solution. Depending on the parametrization, each parameter has a well-defined physical meaning. The obtained ODE system can be used for the estimation and control purposes. The model reduction is performed by minimizing the divergence of flows between the original and reduced systems, and we show that this is equivalent to the minimization of the Wasserstein distance derivative. The method is then tested on the heat equation and on the LWR (Lighthill-Whitham-Richards) model for vehicle traffic.


INTRODUCTION
The problems of controlling high-dimensional systems appear naturally in many areas, including urban traffic, flow dynamics, chemical processes, heating processes, etc. Mathematical models of such systems can be described in various ways, for example, partial differential equations, conservation laws, or networks. However, regardless of the way it is modelled, the problem of controlling the state of the entire system is usually highly complex. It is often preferable to control not the entire state of the system, but some aggregated characteristics, see Klickstein et al. (2016); Casadei et al. (2018). Thus, it is possible to reduce the amount of energy required for control and to reduce sensing requirements.
A typical example of controlling aggregated characteristics is controlling the average state of the system, as in Nikitin et al. (2019). When it is necessary that all system states are close to the average state, variance minimization methods can also be applied (Andreasson et al. (2017); Nikitin et al. (2020)). However, general methods for controlling higher moments of the system require solving the problem of moments closure (Kuehn (2016)). In a particular case, this problem is solved if the system is homogeneous, that is, if the evolution equation for each state equally depends on other states (Zhang et. al. (2018)). For homogeneous systems, mean field analysis (Bick et. al. (2019)) and population density methods (Omurtag et. al. (2000)) can be used.
Other aggregated characteristics of the state of the system may be parameters that describe the spatial properties of the solution. For example, if there is a clear peak in the solution, it would be desirable to be able to describe the dynamics of the position and size of the peak. Often when describing the state of the system it is enough to know a simplified shape of the solution described by several parameters. Moreover, depending on various tasks, various basic shapes may be assumed.
In this work, we present a model reduction method for conservation laws, where the model is reduced to the dynamics of user-defined aggregated characteristics that describe the simplified shape of the solution. Conservation laws are an important class of systems as they can describe various real processes. For example, road traffic is often modelled by a hyperbolic conservation law called LWR model (Lighthill and Whitham (1955); Richards (1956)), and heat distribution is modelled by the parabolic heat equation.
For the model reduction of the parabolic PDEs, the Galerkin approximation (see Li and Qi (2010)) is often used. In this method, the equation is projected onto a set of basis functions, then a finite subset of these functions is selected, and then the final ODE system for projection gains is constructed. For a recent work on controlling PDEs using the Galerkin method and B-splines see Tol et. al. (2019). The Galerkin method is applicable also for nonlinear systems, however, the process of model reduction is linear. The state vector of the obtained ODE system in the general case does not have an explicit physical meaning, and its dimension often turns out to be very large to describe the solution. Many methods have been proposed to refine the solution and find good basis functions, see for example Baker and Christofides (2000) and Barrault et. al. (2004). Hyperbolic conservation laws can create shocks and discontinuities in a finite time, so conventional projective methods do not work. For their approximation, discontinuous Galerkin methods (Cockburn et. al. (2012)) were developed, which can be easily parallelized for the efficient computation. Nevertheless, the dimension of the state vector in this case is enormous.
In this paper, we propose a novel nonlinear model reduction method, in which just one function is used instead of a set of basis functions. This function describes the form of the solution depending on several parameters. The dynamics of the system turns into the dynamics of the shape parameters. The resulting system can be used for estimation and control tasks. We have also shown that the model reduction process minimizes the derivative of the Wasserstein distance (Villani (2008)) between the original and reduced systems.
In Section 2 the general derivation of the reduced model is presented. Section 3 presents one relevant choice of parametrization of the solution shape. Then, an example of applying our method to the reduction of the LWR model is shown in Section 4. Section 5 suggests a method for approximating boundary conditions and gives more examples. Finally, Appendix A is devoted to the analysis of the Wasserstein distance between the original and reduced systems.

PROBLEM FORMULATION AND FORMAL SOLUTION
Let the original system be described by the following model: where the state is the density ρ(x,t). The flow is described by φ (x,t) = φ (x, ρ(x,t), ρ x (x,t)) and can depend on the position, density or its derivative.
We aim to create a reduced system, which is also a conservation law: whereρ(x,t) is the approximated density andφ (x,t) is the approximated flow.
At each time we set the approximated density to have the form ρ(x,t) = g(x, θ (t)), where g(x, θ ) is a function which describes the desired shape based on m parameters θ ∈ R m . This function is assumed to be continuous both in x and θ . As an example one can imagine g(x, θ ) being a Gaussian kernel with θ = (µ, σ ), where µ is the position of the peak and σ is the standard deviation. The parameters θ will constitute the state of the reduced system.
Therefore, we assume that there exists some ODE system which drives the dynamics of the parameters θ : (3) Fromρ(x,t) = g(x, θ (t)) it is possible to write a time evolution equation for the approximated density: By the continuity of g(x, θ ) the time derivative of the approximated density is also continuous. We can imagine that change of the density was caused by some flow which we call an approximated flow. To satisfy (2), this flow should obey a conservation law: Taking an integral and substituting (4), we obtain whereφ 0 (t) is an integration constant and does not affect the dynamics. Finally, assume at initial time point t 0 the density ρ(x,t 0 ) of the original system is equal to the density of the reduced systemρ(x,t 0 ). Now we are ready to define the model reduction procedure. If both conservation laws start from the same initial condition, the natural way to minimize the difference between them is to minimize the L p -difference between the flows, which is shown in the Appendix A to be related to the Wassershein distance derivative minimization in case p = 1. However, for computational purposes we prefer to choose p = 2, which leads to the least squares minimization. Namely, the dynamics F(θ ) for the reduced system can be found as where with B ⊂ R being the domain of the system (or a region of interest).
The flow φ (x,t) in general depends on ρ(x,t), but by our assumption at time point t 0 the density of the original system is the same as the approximated density, thus vector ξ = ( f T ,φ 0 ) T , and write With this notation the cost functional can be written as and its minimization is performed by setting ∂ J ∂ ξ = 0. The minimization of the quadratic function is achieved by solving the linear equation H(θ )ξ = ψ(θ ), and the solution is just Finally, the optimal dynamics for the reduced system iṡ Note that knowing the flow φ (x, ρ, ρ x ) and class of functions g(x, θ ) one can compute H(θ ) and ψ(θ ) symbolically, thus obtaining a closed-form solution to the problem. Moreover, the matrix H(θ ) depends only on the parametrization g(x, θ ) and not on the particular flow φ (x, ρ, ρ x ), therefore it is necessary to symbolically compute it (and its inverse) only once for each chosen parametrization.

A CLASS OF PARAMETRIZATION FUNCTIONS
The most important question which arises while designing the reduced system is the choice of the class of reduced solutions g(x, θ ). One possible solution which is known as Galerkin projection is to take a countable set of basis functions and to define θ to be their multipliers. The examples are the set of all polynomials or the set of all harmonic functions. However this leads to a large number of parameters θ which need to be maintained, especially when the density profile cannot be easily described as a finite sum of basis functions.
What do we suggest is to find a parametrization for each particular case. Here we will present a simple example of one class of functions which can be useful when describing some density with single peak or spike. Such functions are difficult to describe using traditional basis functions like sinusoids or polynomials.
Let θ = (γ, µ, k 1 , k 2 , c 1 , c 2 ), where the meaning of the parameters is as follows: γ is the height of the peak, µ is the position of the peak, k 1 is the slope to the left, k 2 is the slope to the right, c 1 is the constant level to the left, c 2 is the constant level to the right. We can define a piecewise-linear function g(x, θ ) as follows: where p 1 = µ − (γ − c 1 )/k 1 and p 2 = µ − (γ − c 2 )/k 2 . This parametrization is depicted in Fig. 1.
Note that contrary to the Galerkin projections such parametrization can lead to the situation when det H(θ ) = 0, therefore the system (14) can no longer be solved. This happens for example when c 1 = γ. It is clear that in this case the shape becomes degenerate, the parameters become dependent and it is no longer possible to resolve which one should be varied to give the smallest flow discrepancy. Thus the system works as long as it preserves its shape, which is rather obvious if one thinks that the particular class of functions was chosen based on the assumed shape of the real density.

APPLICATION TO LWR SYSTEM
We will show the capabilities of our method on the example of the LWR system: where φ (ρ(x,t)) = ρ(x,t)(2 − ρ(x,t)). This choice of the flow corresponds to the Greenshilds fundamental diagram (Greenshields et. al. (1935)) with ρ max = 2 and v max = 2. This system models the flow of cars on a highway in assumption that the velocity of each car decreases linearly with the density of vehicles nearby.
The results of the comparison are shown in Fig. 2. The original system (16) was numerically solved using Godunov method (Godunov (1959)). Matrix H(θ ) and vector ψ(θ ) were symbolically computed using MATLAB Symbolic Toolbox, and the system (14) was numerically solved using Euler method. We used a space grid with 200 cells for the original system simulation, and the number of time steps was 400 for both systems. We calculated the time needed to simulate both systems: on average, simulation of the original system took 0.207241 seconds, and simulation of the reduced system took 0.027709 seconds, thus being almost 10 times faster. We believe that by designing a specialized software instead of using a general Toolbox one can achieve much higher performance. Figure 2. Comparison of the real and reduced solutions to the LWR equation (16). Insets a) and b) show the real density ρ(x,t) and approximated densityρ(x,t) respectively. Inset c) shows the initial density ρ(x, 0) with black dashed line, the real density ρ(x,t) at t = 2 with blue line and the approximated densityρ(x,t) at t = 2 with red line. Inset d) shows the same data for t = 8.
It is clear that the reduced solution perfectly tracks the position and the slopes of the peak, and the difference between the real density and the approximated density arises only because of the non-smoothness of the reduced solution.
It is also interesting to note that if one continues simulation further, the reduced system will fail at time moment t = 10.1, because the shock arises on the left slope and k 1 becomes infinity. It is a known property of the LWR system which can produce shocks in a finite time. In our assumption the function g(x, θ ) should be continuous both in x and θ for the correct definition of the artificial flow (6).

BOUNDARY PROBLEMS
Up to now not a word was said about boundary conditions which affect the solution to the original system and which should be taken into account properly in the reduced system.
Essentially all the analysis performed in the previous sections was based on the assumption that the solutions evolve in R, with only exception being the cost functional (8). Therefore boundary conditions were not taken into account either in original or reduced systems. Now assume that the original system is given by the equation (1) defined on the domain B = (a, b) ⊂ R for some a < b. Let one of the boundary flows φ (a,t) = φ in (t) or φ (b,t) = φ out (t) (or possibly both of them) be given. The flows can be either given explicitly, as for example in the case of heat generators in the heat equation, or they can depend on the density itself, as in the case of demand and supply flows in the LWR system.
The given boundary flows work as constraints for the flow discrepancy minimization problem. Namely, if the inflow φ in (t) is given, the solution ξ to the problem of minimization of the cost functional J should satisfy the constraint and similarly for outflow in case φ out (t) is given. Figure 3. Comparison of the real and reduced solutions to the heat equation (22) with φ in = 0 and φ out = 0. Insets a) and b) show the real density ρ(x,t) and approximated densityρ(x,t) respectively. Inset c) shows the initial density ρ(x, 0) with black dashed line, the real density ρ(x,t) at t = 2 with blue line and the approximated densityρ(x,t) at t = 2 with red line. Inset d) shows the same data for t = 8.
The constraints can be written in a unified manner if one defines matrix C(θ ) and column-vector d such that they have one row if only one condition is given and two rows if both boundary conditions are set. In the latter case Now the constrained minimization of J is equivalent to the minimization of the Lagrangian and the solution to this problem, omitting the dependency on θ , is given by (20) From ξ the dynamics of θ can be easily recovered by discarding the last row:θ = H −1 ψ −C T λ 1...m .
(21) By putting the constraints on both boundaries one guaranties the conservation of mass, therefore the overall mass in the original and the reduced system are always equal.

Applications to heat equation and LWR model
As an example one can look at the heat equation: where φ (ρ(x,t)) = −ρ x (x,t), and we consider its approximation using (15). The simulations were performed with the same space and time discretization as in the previous section.  (22) with φ in = 0.02 and φ out = −0.08. Insets a) and b) show the real density ρ(x,t) and approximated densityρ(x,t) respectively. Inset c) shows the initial density ρ(x, 0) with black dashed line, the real density ρ(x,t) at t = 2 with blue line and the approximated densityρ(x,t) at t = 2 with red line. Inset d) shows the same data for t = 8.
Assuming that there is no heat transfer between the domain and the environment, we set the boundary conditions φ in = φ out = 0. The results are presented in Fig. 3, and it is seen that the reduced model approximates the original one almost perfectly.
Alternatively, we can set non-zero boundary conditions for the heat equation, for example φ in = 0.02 and φ out = −0.08, which corresponds to the injection of heat into the system from both boundaries. The results are shown in Fig. 4. It is clear that although the reduced system cannot capture the shape of the real solution near the right boundary at time t = 8, it still does its best in terms of the chosen shape.
The control over one boundary can be demonstrated on the LWR example. Assume we set the inflow φ in (t) = 0.3 + 0.15 sin(3t). Then the reduced system will average the high frequency components, while tracking the initial peak. The results are presented in Fig. 5. Figure 5. Comparison of the real and reduced solutions to the LWR equation (16) with φ in (t) = 0.3 + 0.15 sin(3t). Insets a) and b) show the real density ρ(x,t) and approximated densityρ(x,t) respectively. Inset c) shows the initial density ρ(x, 0) with black dashed line, the real density ρ(x,t) at t = 2 with blue line and the approximated densityρ(x,t) at t = 2 with red line. Inset d) shows the same data for t = 8.

CONCLUSION
In this paper we presented a method of describing 1D conservation law system based on a notion of a solution shape. We reduced system state to a set of well-tractable shape parameters and derived their dynamics, providing closed-form solution.
The reduced model can potentially be used for the control design based on the aggregated characteristics of the system.
Throughout this paper we assumed that the initial condition for the original system belongs to the class g(x, θ ) and thus the approximation error at the beginning is zero. For the future research it remains to develop the methods for the initial state projection and also to determine bounds on the solution approximation error in a long run.
Further, to fight with degeneracy of the solutions when the shape is violated it is interesting to develop a methods for reparametrizations and changes of shape, as well as to allow discontinuous shape functions. It is possible to show that the minimization of flow discrepancy leads to the minimization of Wasserstein distance divergence between real and reduced solution.
The L p -Wasserstein distance between two nonnegative densities ρ 0 (x) and ρ 1 (x) of equal mass is defined as where T is the set of all possible transformations over the domain B that transfer the mass from one configuration to another. Namely, We will show that the L 1 -minimization of flows is equivalent to the minimization of the time derivative of L 1 -Wasserstein distance.
Assume at some time moment t 0 the state of the original system ρ(x,t 0 ) and the reconstructed state of the reduced system ρ(x,t 0 ) are equal. Then the L 1 -Wasserstein distance is zero and the transformation T which achieves minimum is identity, T (x) = x. Equivalently we can define time-dependent transformation T (x,t) such that T (x,t 0 ) ≡ id. Now take the time derivative of the L 1 -Wasserstein distance for this particular transformation: where we used the fact that ρ(x,t 0 ) 0.