Inverse demand tracking in transportation networks

Göttlich, Simone; Mehlitz, Patrick; Schillinger, Thomas Article — Published Version Inverse demand tracking in transportation networks Mathematical Methods of Operations Research Provided in Cooperation with: Springer Nature Suggested Citation: Göttlich, Simone; Mehlitz, Patrick; Schillinger, Thomas (2024) : Inverse demand tracking in transportation networks, Mathematical Methods of Operations Research, ISSN 1432-5217, Springer, Berlin, Heidelberg, Vol. 100, Iss. 3, pp. 635-668, https://doi.org/10.1007/s00186-024-00875-y This Version is available at: https://hdl.handle.net/10419/314974 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ Mathematical Methods of Operations Research (2024) 100:635–668 https://doi.org/10.1007/s00186-024-00875-y ORIGINAL ARTICLE Inverse demand tracking in transportation networks Simone Göttlich1·Patrick Mehlitz2·Thomas Schillinger1 Received: 14 July 2023 / Revised: 9 August 2024 / Accepted: 15 August 2024 / Published online: 5 October 2024 © The Author(s) 2024 Abstract This paper deals with the reconstruction of the desired demand in an optimal control problem, stated over a tree-shaped transportation network which is governed by a linear hyperbolic conservation law. As desired demands typically undergo fluctuations due to seasonality or unexpected events making short-term adjustments necessary, such an approach can exemplary be used for forecasting from past data. We suggest to model this problem as a so-called inverse optimal control problem, i.e., a hierarchical optimization problem whose inner problem is the optimal control problem and whose outer problem is the reconstruction problem. In order to guarantee the existence of solutions in the function space framework, the hyperbolic conservation law is interpreted in weak sense allowing for control functions in Lebesgue spaces. For the computational treatment of the model, we transfer the hierarchical problem into a nonsmooth single-level one by plugging the uniquely determined solution of the inner optimal control problem into the outer reconstruction problem before applying techniques from nonsmooth optimization. Some numerical experiments are presented to visualize various features of the model including different types of noise in the demand and strategies of how to observe the network in order to obtain good reconstructions of the desired demand. Keywords Inverse optimal control ·Linear hyperbolic conservation laws · Transportation networks Mathematics Subject Classification 49J20 ·65M32 ·90C33 ·90C35 BSimone Göttlich [email protected] Patrick Mehlitz [email protected] Thomas Schillinger [email protected] 1School of Business Informatics and Mathematics, University of Mannheim, 68159 Mannheim, Germany 2Department of Mathematics and Computer Science, Philipps-Universität Marburg, 35032 Marburg, Germany 123 636 S. Göttlich et al. 1 Introduction Flow problems over energy and supply networks model a broad range of interesting applications, see (Bressan et al. 2014) for a survey. In this paper, we investigate transportation networks of tree shape where the flow on edges is modeled, for simplicity, via (linear) hyperbolic conservation laws, as typically used for electric transmission lines (Göttlich et al. 2016), heating networks (Rein et al. 2020), or networks of gas pipelines (Banda et al. 2006; Gugat et al. 2018). A control function is used to model the inflow at some source vertex, and the aim of optimization is to choose this function in such a way that certain desired demands at the sinks of the network are tracked as close as possible. As mentioned in some recent contributions, see (Göttlich et al. 2019; Göttlich and Schillinger 2022a,b), these desirable demands are subject to perturbations, noise, or other sources of stochasticity. In the aforementioned papers, this issue has been faced by modeling the problem as a stochastic optimal control problem which is influenced by randomness via appropriately chosen stochastic processes. In this paper, we are concerned with related phenomena. Let us consider the following practically relevant situation. There exists a company (C2) which appoints a second company (C1) to deliver a certain amount of electricity/heat/gas at the demand vertices over time by inserting the requested product at the source of the network over time. In this regard, C1 has to solve the optimal control problem mentioned above. We now enrich the considered situation by assuming that there is a network operator (NO), different from C1 and C2, which partially observes the flow along the network and, depending on this, charges C1 and C2 to pay some tax for employing the network. As outlined above, the desired demands requested by C2 are subject to stochastic influences and, additionally, may vary due to a seasonal behavior. From past data, NO now wants to forecast the desired demand of C2 and the associated actions of C1, exemplary for fixing taxes to plan future income. Typically, NO is not aware of the desired demand as he only observes the actual network flow along some but, most likely, not all edges of the network (as it might be expensive to equip the overall network with sensors or to run them on each edge over all time). Furthermore, the forecasting model should be capable of recognizing seasonal behavior of the desired demands as it is exemplarily presented for an electricity market in Coskun and Korn (2021). In order to model this situation, we consider it from the viewpoint of inverse optimal control, i.e., we aim to identify parameters in an optimal control problem (and not only in a dynamical system). Here, the optimal control problem of interest is the aforementioned network flow problem, and the appearing desired demand plays the role of this parameter. We assume that we are given observed (but, most likely, noisy) pairs of optimal inflow and optimal network flow, and aim to reconstruct the desired demands which are modeled as a convex combination of given ansatz functions. It is, thus, our goal to find the associated weight parameters which characterize a suitable standard (periodically emerging) choice for the desired demand. As we are interested in the robustness of our approach, we consider additional perturbations in the model and study different types of temporal restrictions in the observation of the network to evaluate whether these are sufficient for good forecasting. Naturally, the model of interest is a hierarchical optimization problem with two decision levels. Coming back to our exemplary situation from above, at the outer 123 Inverse demand tracking in transportation networks 637 (or upper-level) problem, the NO is in position to partially observe the network and chooses certain weights, which then give a tangible desired demand. At the inner (or lower-level) problem, C1 now can solve the network flow problem. Along those parts of the network, which are observed by NO, the latter can compare the past data and the real-time data obtained from the inner problem for this particular choice of the weight parameters. Noting that this decision order leads to a well-posed problem, NO aims to choose the weight parameters in such a way that past data and real-time data match as good as possible. As our model has two decision levels, it is a so-called bilevel optimization problem. For more than 50 years, bilevel optimization is a major field of research in mathematical programming due to numerous underlying applications e.g. in data science, economy, finance, machine learning, or natural sciences, see (Bard 1998; Dempe 2002; Shimizu et al. 1997) for an introduction and Dempe (2020) for a recent survey which presents an overview of contributions in this area. Recently, bilevel optimization turned out to be of particular interest in the context of transportation or energy networks, see e.g. Dempe et al. (2015). This also includes the rapidly growing field of hierarchical control, see e.g. Mehlitz and Wachsmuth (2020) for an overview, and, particularly, so-called inverse optimal control already mentioned earlier, see (Hinze et al. 2009; Tröltzsch 2010; Troutman 1996; Vinter 2010) for an introduction to the topic of optimal control. Inverse control possesses several interesting applications e.g. in the context of human locomotion, see (Albrecht et al. 2012; Albrecht and Ulbrich 2017; Albrecht et al. 2010; Mombaur et al. 2010). The theory on inverse optimal control including ordinary and partial differential equations addresses the existence of solutions, optimality conditions, and solution algorithms, see e.g. (Dempe et al. 2019; Friedemann et al. 2023; Harder and Wachsmuth 2019;Hatzetal.2012; Holler et al. 2018; Suryan et al. 2016) and is developing fast. In abstract bilevel optimization, two decision makers, a leader and a follower, need to choose variables in order to minimize their associated cost function which also depends on the variables of the other decision maker, respectively. More precisely, the leader chooses his variables first which are handed over to the follower who now can solve his optimization problem (which is parametric in the leader’s variable) to global optimality. The solutions are then given to the leader, who now can evaluate his objective. Often, one assumes that leader and follower cooperate in order to optimize the leader’s objective, and this procedure is referred to as the optimistic approach to the problem, see (Zemkoho 2016) for an overview of other approaches avoiding ill-posedness in bilevel optimization. The leader’s and follower’s problem are often referred to as upperand lower-level problem, respectively. As the follower has to determine globally optimal solutions of his problem by nature of bilevel optimization, one typically requires that the lowerlevel problem is convex in the follower’s variable in order to circumvent issues related to nonconvex global optimization at the lower-level stage. We start our investigations by modeling the problem of interest as an inverse control problem in Sect.2. Therefore, we first study the existence of solutions for linear hyperbolic conservation laws in a function space which is suitable for optimal control before setting up the lowerand upper-level problem consecutively. Furthermore, we demonstrate that the resulting optimization problem possesses an optimal solution in the function space setting we are investigating. In Sect. 3, we address the computa123 638 S. Göttlich et al. tional treatment of the model. Section3.1 describes our approach to the numerical solution of the problem. As it is analytically possible to compute the network flow associated with the input, we are in position to distill a state-reduced version of the parametric optimal control problem. The associated solution operator, which, at least in pointwise fashion, can similarly be computed analytically due to the nice structure of the problem, turns out to be a nonsmooth single-valued mapping. Plugging the latter into the superordinate reconstruction problem and performing a suitable discretization, we end up with a nonsmooth optimization which we solve with the aid of MATLAB’s patternsearch solver in default mode. The general set-up of our computational experiments is carved out in Sect. 3.2. Numerical results are presented in Sect.3.3 in order to visualize the effectiveness and several different features of the approach. Particular focus is laid on the robustness of the model with respect to additional uncertainties, restricted observation options, and the presence of additional inflow constraints. Some concluding remarks close the paper in Sect. 4. 2 The model problem In this section, we set up the model of our interest. First, we discuss the particular shape of the lower-level parametric optimal control problem in Sect.2.1. Therefore, we first present the underlying network dynamics and discuss regularity features of associated solutions. Second, the lower-level objective function is constructed, and solvability of the overall lower-level problem is discussed. In Sect. 2.2, we derive the superordinate upper-level problem and demonstrate that it possesses an optimal solution in the function space setting. 2.1 The lower-level problem In this subsection, we are concerned with the derivation and analysis of the lowerlevel optimal control problem. To start, we state the lower-level dynamics and discuss existence and uniqueness of solutions associated with this system. Afterwards, we set up the (parametric) lower-level problem, show that, for each set of parameters, it possesses a unique solution, and investigate properties of the associated solution operator. 2.1.1 Setting up the network and network dynamics We consider a directed graph G=(V,E)which is a tree (in the sense that whenever the directed edges are interpreted as undirected, then the resulting graph would be free of cycles). Let us use the notation V:= {v0,...,v n}and note that |E|=nby nature of trees. Some more details on Gand the notation we are going to exploit are discussed below. •The uniquely determined source vertex of the network Gis v0∈V. Furthermore, we assume that v0is a leaf of G, i.e., there is only one edge which leaves v0, and the vertex at its end will be denoted by v1. 123 Inverse demand tracking in transportation networks 639 Fig. 1 An exemplary network with VD={v4,v 5,v 6},VI={v1,v 2,v 3},E+(2)={(4), (5)},and ED={(4), (5), (6)} •In VD⊂V, we collect all vertices which possess no outgoing edges. These are the demand vertices. •All remaining intermediate (or inner) vertices of the network are collected in the set VI:= V\(VD∪{v0}). •For vi∈V\{v0}, we identify the uniquely determined edge which ends at viby (i). •The set E+(i)is used to denote the set of all edges starting at vertex vi. Furthermore, we use ED:= {(i)∈E|vi∈VD}to denote the set of edges that end at a demand vertex. Clearly, |ED|=|VD|. We visualize the above notation in Fig.1. For the theory of this paper, it is not mandatory that the vertex v0possesses just one outgoing edge. One can interpret v0as an upstream supersource. Besides, this additional assumption simplifies the notation because we can abstain from the introduction of distribution parameters at the inflow vertex later on. At the source v0, the injection of flow over time T:= (0,T), where T>0isthe final time, is modeled by the control variable u:T→Rwhich has to be chosen from an appropriate function space. The flow over (i)at time t∈Tat the spatial coordinate x∈will be denoted by z(i)(t,x). Here, we assume that := (0,ω)is a bounded real interval. The density has to obey the linear hyperbolic conservation law z(i) t(t,x)+λ(i)z(i) x(t,x)=0,a.e. on T×, (i)∈E,(2.1a) z(i)(0,x)=0,a.e. on , (i)∈E,(2.1b) λ(1)z(1)(t,0)=u(t), a.e. on T,(2.1c) λ(k)z(k)(t,0)=αi,kλ(i)z(i)(t,ω), a.e. on T,v i∈VI,(k)∈E+(i). (2.1d) Particularly, the flux functions of the conservation law are of linear structure. For each i∈{1,...,n},λ(i)>0 is a given constant. Above, for each vi∈VIand (k)∈E+(i), αi,k>0 is a constant such that (k)∈E+(i)αi,k=1 holds, i.e., the coefficients αi,k model how the flow splits at vertex viinto the flows along the edges from E+(i). This way, (2.1d) conserves the flow. We note that the theory can be extended to more general situations. Exemplary, standard linear damping terms of type μ(i)z(i)(t,x) 123 640 S. Göttlich et al. can be incorporated in (2.1a) for real constants μ(i)>0 for each (i)∈Ewithout any problem. Under additional assumptions, the coefficients λ(i)and μ(i)may also depend on time. Without loss of generality one could choose ω:= 1. However, in order to clearly distinguish between temporal and spatial variables in notation, we stick to the seemingly more general situation where ω>0 is arbitrary. Furthermore, the findings in this paper extend to connected networks without cycles, but apart from a more difficult notation, which also allows for vertices where flows are merged, we do not believe that such a model comes along with a significantly different theory. We, thus, concentrate on tree-shaped networks. 2.1.2 Discussion of the hyperbolic conservation law Let us first review a classical existence result for the linear hyperbolic conservation law (2.1). Therefore, we define a suitable control space by C1 00(T):= u∈C1(T)|u(0)=0,u(0)=0. We equip C1 00(T)with the classical C1-norm, and note that this space is a closed subspace of C1(T). The proof of the following result, which is based on the method of characteristics, can be distilled from Bressan (Bressan 2000, Section 3.1, Theorems 3.4 and 3.6) under the condition that we only consider positive velocities on the network and, thus, all waves are moving with positive speed. Proposition 2.1 For each u ∈C1 00(T), the hyperbolic conservation law (2.1) possesses a unique solution z := (z(1),...,z(n))∈C1(T×, Rn). The latter is explicitly given by ∀(t,x)∈T×:z(1)(t,x)=1 λ(1)u(t−x/λ(1))t−x/λ(1)>0, 0t−x/λ(1)≤0(2.2) on edge (1), and for each i ∈{1,...,n}such that vi∈VIand (k)∈E+(i), we find ∀(t,x)∈T×:z(k)(t,x)=αi,kλ(i) λ(k)z(i)(t−x/λ(k),ω) t−x/λ(k)>0, 0t−x/λ(k)≤0. (2.3) Additionally, there is a constant κ>0, not depending on u, such that zC1(T×,Rn)≤ κuC1(T). Let us note that formula (2.3) can be used recursively to determine the solution along all edges of the network. Indeed, based on (2.2), the solution along all arcs from E+(1)can be computed. Next, using (2.3), it is possible to determine the flow along all edges starting in those vertices which are the end vertex of some edge in E+(1). Repeating this procedure, one can iterate through the whole network. 123 Inverse demand tracking in transportation networks 641 Clearly, Proposition 2.1 justifies to introduce a map from C1 00(T)to C1(T×, Rn) which assigns to each control function from C1 00(T)the associated uniquely determined solution of (2.1). This mapping is a linear operator which is continuous by Proposition 2.1. Since we are interested in the optimal control of the system (2.1), working with the control space C1 00(T)induces some inherent difficulties. First, this space is nonreflexive, i.e., to show the existence of optimal solutions for optimization problems over (2.1) and the superordinate inverse optimal control problem, which we state in Sect.2.2, would be challenging. Second, the dual of this space, which naturally arises when using the adjoint approach for the derivation of optimality conditions, is large and difficult to handle numerically. It is, thus, a reasonable task to reconsider (2.1) from the viewpoint of control functions u∈L2(T). Besides, this choice allows for discontinuous controls which can be exploited to model switches in the inflow. Observe that (2.1) does not need to possess a classical solution in the sense of Proposition 2.1 anymore whenever the control function is not continuously differentiable. To proceed, we follow (Keimer 2014, Section 2.2), see (Gugat et al. 2015, Section 2) as well, to introduce a suitable weak formulation of (2.1) as stated below. First, for the state z(1), we demand Tτ z(1)(t,x)(ϕt(t,x)+λ(1)ϕx(t,x))dxdt =−Tτ u(t)ϕ(t,0)dt∀ϕ∈Wτ(2.4) for all τ∈T, where Tτ:= (0,τ)and Wτ:= ϕ∈C1(Tτ×)  ϕ(·,ω)=0on Tτ ϕ(τ,·)=0on is the space of test functions. Similarly as above, we demand Tτ z(k)(t,x)(ϕt(t,x)+λ(k)ϕx(t,x))dxdt =−αi,kλ(i)Tτ z(i)(t,ω)ϕ(t,0)dt∀ϕ∈Wτ(2.5) for all τ∈T,vi∈VI, and (k)∈E+(i). A function z∈C(, L2(T,Rn)) satisfying these requirements is referred to as a weak solution of the hyperbolic conservation law (2.1). Recall that the function space C(, L2(T,Rn)) comprises all functions z:T×→Rnsuch that, for each x∈,z(·,x)belongs to L2(T,Rn), and x→ z(·,x)∈L2(T,Rn)is continuous. Let us emphasize that the boundary conditions (2.1b), (2.1c) are incorporated in this alternative formulation of the dynamics also in weak sense only (by definition of the space Wτ) since pointwise considerations are meaningless in Lebesgue spaces. 123 642 S. Göttlich et al. The following result shows that the (classical) solution characterized in Proposition 2.1 (with controls chosen from C1 00(T)) also provides the uniquely determined weak solution of the hyperbolic conservation law (2.1) if the control is chosen from L2(T). Proposition 2.2 For each u ∈L2(T), the function z := (z(1),...,z(n))∈ C(, L2(T,Rn)) characterized via (2.2), (2.3) is the uniquely determined weak solution of the hyperbolic conservation law (2.1). Additionally, there is a constant κ>0, not depending on u, such that zC(,L2(T,Rn)) ≤κuL2(T). Proof Let us start to show that z(1)given in (2.2) satisfies (2.4) for each τ∈Tand given u∈L2(T). Therefore, we introduce a function ¯u∈L2((−ω/λ(1),T)) by ∀t∈(−ω/λ(1),T):¯u(t):= u(t)t>0, 0t≤0. Using a coordinate transformation with respect to the new domain τ:= {(s,x)∈R2|x∈, s∈(−x/λ(1),τ −x/λ(1))}, we find, for each ϕ∈Wτand ¯ϕ(s,x):= ϕ(s+x/λ(1),x)for all (s,x)∈τ,the identities Tτ z(1)(t,x)(ϕt(t,x)+λ(1)ϕx(t,x))dxdt =1 λ(1)Tτ ¯u(t−x/λ(1))(ϕt(t,x)+λ(1)ϕx(t,x))dxdt =τ ¯u(s)¯ϕx(s,x)d(s,x) =τ −ω/λ(1)min(ω,λ(1)(τ−s)) max(0,−λ(1)s) ¯u(s)¯ϕx(s,x)dxds =τ 0 u(s)min(ω,λ(1)(τ−s)) 0 ¯ϕx(s,x)dxds =τ 0 u(s)( ¯ϕ(s,min(ω, λ(1)(τ −s))) −¯ϕ(s,0))ds =−Tτ u(s)ϕ(s,0)ds. Above, we used the fact that the determinant of the Jacobian associated with the chosen coordinate transform is 1, the fundamental theorem of calculus, and ¯ϕ(s,min(ω, λ(1)(τ −s))) =ϕ(τ,λ(1)(τ −s)) =0λ(1)(τ −s)<ω, ϕ(s+ω/λ(1),ω)=0λ(1)(τ −s)≥ω, 123 Inverse demand tracking in transportation networks 649 Let us demonstrate that the reconstruction problem (UL) possesses an optimal solution. Proposition 2.9 The optimization problem (UL) possesses a globally optimal solution. Proof We note that (UL) can be transferred into a finite-dimensional optimization problem by plugging the lower-level solution operator into the objective function. It is obvious that a point β∈(Rm)|VD|is a global minimizer of the resulting controlreduced problem if and only if (β, (β)) is a global minimizer of (UL). By continuity of , see Proposition 2.7, and continuity of Cas well as D, the objective function of the reduced problem is then continuous, while its feasible set (m)|VD|is nonempty and compact. Thus, the reduced problem possesses a global minimizer ¯ β∈(Rm)|VD|by the Weierstraß theorem, and this yields that (¯ β,( ¯ β))solves (UL) to global optimality.  Although being globally Lipschitz continuous, see Proposition 2.7, the lower-level solution operator , which, at least in discretized form, see Appendix A, can be represented as the composition of a linear, continuous operator and the projection onto Uad, is likely to be nonsmooth apart from the special situation where no control constraints are present, see Remark 2.8. Eliminating the control variable uin (UL) by plugging into the objective function, thus, leads to a finite-dimensional but nonconvex, nonsmooth optimization problem with polyhedral constraints. Whenever Uad =L2(T)holds, is linear, see Remark 2.8 again, and (UL) is actually a convex optimization problem. In this particular situation, numerical methods which identify stationary points of (UL) may already compute global minimizers of the problem. This is a rare property in hierarchical optimization where the multilevel structure is, typically, a source of nonconvexity and nonsmoothness, and this problem we also face in the general setting where control constraints are present. 3 Numerical solution and computational results In this section, we first describe how (UL) can be solved in numerical practice. Second, results of some computational experiments are presented. 3.1 Numerical solution of the problem For the network discretization, we choose a time grid (tj)J j=1of J∈Ndiscretization points such that tj:= (j−1)tfor all j∈{1,...,J}, where t>0isagiven temporal stepsize, and a spatial discretization of each edge (i), represented by the interval (0,ω),as(x(i) q)L(i) q=1, where L(i)∈Nis the number of discretization points, x(i) q:= (q−1)x(i)for all q∈{1,...,L(i)}, and x(i)>0 is the spatial stepsize for edge (i). The transported quantities z(i) j,qat time tjand position x(i) qgiven by the PDE in (2.1a) are calculated using a left-sided upwind scheme, i.e., z(i) j,q=z(i) j−1,q−t x(i)λ(i)z(i) j−1,q−z(i) j−1,q−1,j∈{2,...,J},q∈{2,...,L(i)}. 123 650 S. Göttlich et al. We also note that (2.1b) translates into z(i) 1,q=0 for all q∈{1,...,L(i)}.Atthe junctions, according to (2.1c) and (2.1d), we require z(1) j,1=uj λ(1),z(k) j,1=αi,k λ(i) λ(k)z(i) j,L(i),v i∈VI,(k)∈E+(i), j∈{1,...,J} where uj:= u(tj)for all j∈{1,...,J}.For t x(i)λ(i)=1, the upwind scheme shows no diffusion. Therefore, we set x(i):= λ(i)twhich leads to different spatial grids on the different edges whenever the respective coefficients λ(i)are not the same. We use this discretization for a finite differences approximation of the lower-level problem (LL(β)). We define S(i),L(i)∈RJ×Jto be the (discrete) realization of S(i) ω such that J ν=1S(i),L(i) j,ν uνapproximates the influence of the discretized inflow on the density z(i) j,L(i)at time tjand spatial point ω. Further, we denote the discrete versions of the demand profiles D(i) 1,...,D(i) mfor edge (i)∈EDby ˜ D(i) 1,..., ˜ D(i) m∈RJ. For our computations, we will exploit that the columns of S(i),L(i)are orthogonal to each other. This is the case since, due to the special structure of the PDEs, there is a one-to-one correspondence between the inflow into the system and the outflow out of the system. Therefore, in the discretized setting, there is a unique time point for the inflow that determines the outflow at the corresponding outflow time. This property enforces the matrix S(i),L(i)to be nonzero on its subdiagonal. Consequently, S(i),L(i) is orthogonal. For a given convex combination of base demands by the vector βand using (2.10), we obtain the optimal inflow in the discretized setting in the absence of control constraints when solving the linear system Au −Bβ=0 where Ais given by A:=  (i)∈EDS(i),L(i)S(i),L(i)+σIJ, where IJ∈RJ×Jis the identity matrix, and B:= Q(i)(i)∈ED. Above, for each (i)∈ED,Q(i)∈RJ×mis given by Q(i):= S(i),L(i)˜ D(i) 1... S(i),L(i)˜ D(i) m. We note that the discretized lower-level problem is equivalent to min u{1 2uAu −(Bβ)u|ua≤u≤ub}, where ua,j:= ua(tj)and ub,j:= ub(tj)for all j=1,...,J. We obtain the solution of this problem by projecting the solution of the linear equation Au −Bβ=0 onto 123 Inverse demand tracking in transportation networks 651 Fig. 2 The network considered in Sect.3.3 with VD={v7,v 8,v 9,v 10,v 11},VI={v1,v 2,v 3,v 4,v 5,v 6}, and ED={(7), (8), (9), (10), (11)} the feasible box, since Ais a diagonal positive definite matrix by orthogonality of S(i),L(i),(i)∈ED, see Appendix Afor details. For the upper-level problem (UL), we apply the same discretization technique with different stepsizes, see Sect. 3.3, and consider, if not specified differently, the observation operator Cin which we only observe the densities at the demand vertices from VD, corresponding to the last discretization points of the edges in ED,aswellasat the first discretization point of edge (1), monitoring the inflow at v0. Additionally, D is the zero operator in our experiments. Further details and some numerical examples are explained in Sect.3.3 where it is also described how Cand Dcan be adjusted. Inserting the discretized solution operator of the lower-level problem into the objective function of the discretized upper-level problem results in a nonsmooth optimization problem with affine constraints, and we solve the latter using MATLAB’s patternsearch solver in default mode. We want to emphasize that the performance of this optimization routine heavily depends on the initial point that is handed over to the solver. This, however, is not surprising as the considered nonsmooth problem of interest is nonconvex and, thus, likely to possess several local minimizers and stationary points which are different from its global minimizers. As the model is designed to reconstruct certain reference parameters from noisy data, we initialize patternsearch with a perturbed version of these reference parameters to face this problem. We note that, in the absence of lower-level control constraints, the resulting single-level problem is a simple convex quadratic problem which can be solved, exemplary, with the aid of MATLAB’s quadprog routine, and the aforementioned issues do not occur. 3.2 General set-up of experiments We consider the tree-shaped network presented in Fig.2in which each edge has a length of ω=1. The velocities are chosen identically for all edges, we use λ(i)=10, i=1,...,11. The stepsizes are given by t=1 60 ,x=1 6for the backward calculation and 123 652 S. Göttlich et al. t=1 70 ,x=1 7for the forward calculation, which are chosen differently to avoid an inverse crime, see (Colton and Kress 2013, page 154), for the unconstrained examples, and t=1 20 ,x=1 2(backward calculation), t=1 30 ,x=1 3(forward calculation) when applying constraints to the inflow in order keep reasonable running times. Note that x(i)=xis exploited, i=1,...,11. In both cases, the Courant– Friedrichs–Lewy condition holds true with equality, i.e., t xλ(i)=1, i=1,...,11, to avoid diffusion in the numerical scheme. The distribution parameters are set to α1,2=0.65,α 2,4=0.7,α 4,7=0.5,α 6,10 =0.4, α1,3=0.35,α 2,5=0.3,α 4,8=0.5,α 6,11 =0.6. We consider the evolution of the demand within one week, i.e., T=168 where one time unit represents one hour and assume four underlying base demand levels which are visualized in Fig.3and chosen as •a time constant level of the demand: D1(t)=4, •a daily varying level at which we attain the highest level in the morning: D2(t)=2+sin (π(t−2)/12), •a daily varying level at which we attain the highest level in the afternoon: D3(t)=2+sin (π(t−10)/12), •a level that illustrates the lower demand during the weekend: D4(t)=1[0,120](t). These choices can similarly be found for example for the electricity market in Coskun and Korn (2021) and describe the identified two-peak pattern of demand in the intraday market (D2,D3) as well as the phenomenon referred to as the weekend effect (D4). For the prototypical demand profiles, we make use of D(i) := ˆ d(i)D, i∈{7,...,11},∈{1,2,3,4}, where ˆ d(7)=0.2275,ˆ d(8)=0.2275,ˆ d(9)=0.195,ˆ d(10)=0.14,ˆ d(11)=0.21. This choice proportionally accounts for the different distribution parameters in the network. The historical observations are basically generated using the initial weights (β1,β 2,β 3,β 4)=(0.2,0.15,0.2,0.45). (3.1) In every time step and for every demand vertex, the base demand levels are perturbed by random variables Z(i) 1∼N(0,1), Z(i) 2∼N(0,1/4), Z(i) 3∼N(0,1/4), Z(i) 4∼N(0,1/4), such that the historically desired demands are given by realizations of D(i) d= 4  =1 βˆ d(i)D+Z(i) ,i∈{7,...,11}.(3.2) 123 Inverse demand tracking in transportation networks 653 Fig. 3 Illustration of the four base demand levels The historically observed pairs (zo,uo)are computed as solutions of the associated problem (2.12). 3.3 Documentation of experiments In the following, we investigate different variants of the bilevel optimization problem (UL). The standard version is presented in Sect.3.3.1, and settings with additional perturbations in the historical observations are shown in Sects.3.3.2 and 3.3.3.A time-restricted observation operator Cis investigated in Sect. 3.3.4. All subsections are constructed in a similar way. First, we present exemplary historical demand observations, then we provide a comparison of the inand outflows for the means of the historical observations and the initially chosen βas well as for the reconstructed βin a framework without an inflow constraint, which can also be considered as a framework with a high constraint that does not really affect the inflow. These illustrations are presented for the inflow vertex and the demand vertex v7(the behavior at all other demand vertices is similar). We can verify that, on the one hand, the optimal inflows are calculated correctly and, on the other hand, see whether the reconstruction of the weights βwas successful. The second aspect is further underlined by a table presenting the means and variances for βof a Monte Carlo simulation of N=40 runs for 123 654 S. Göttlich et al. different numbers of historical observations p. Second, we repeat the investigations of each subcase based on a medium inflow constraint ub≡2 and a low inflow constraint ub≡1.5, where we also ensure nonnegative inflows, i.e., ua≡0, the latter being nonrestrictive as the desired demand at the vertices in VDis nonnegative. 3.3.1 Standard model without additional adjustments In this scenario, no further perturbations or model changes are included, and we consider the framework presented in the previous sections. Three examples for historical observations are given in Fig.4which show the sinusoidal behavior of demand, as well asthedropfort>120 during the weekend. Furthermore, we detect the stochastic noise in the demands, however, still verify that the demands show a very similar structure. The comparison of the inflow and outflow for demand vertex v7are presented in Fig.5, where the blue curve shows the mean values of the p=6 historical observations, the yellow dotted line represents the curve for the true βgiven in (3.1), and the red line the inor outflow for the reconstructed β. All considerations were made without constraining the inflow control. It can be concluded that all three curves match very well, which means that, on the one hand, the inflow is calculated appropriately and, on the other hand, also the weights of the base demands are reobtained very well. The outflow behavior at the demand vertices v8,...,v 11 shows similar patterns and is (for brevity of presentation) not illustrated. At the beginning and the end of the considered time horizon, some curves in Fig.5decay to zero or show a jump. This can be explained by the fact that around time t=0, it takes some time until (starting from an empty system) the first inserted quantity reaches the demand vertex. Therefore, the outflows are zero in the very beginning of the time period. Conversely, the inflow for times close to T=168 vanishes, since these quantities do not reach the demand nodes within the considered time horizon. The increase at T=168 in the outflow figure can be explained by considering T=168 to be Monday already, where the demand is larger again. Similar artifacts show up in some other figures in this section due to analogous reasons. Table 1shows the means and variances of the reconstructed weights for the base demands for different numbers of perturbed historical observations in a Monte Carlo simulation of N=40 runs and underlines the results from Fig. 5quantitatively. As it can be expected for larger numbers of historical observations, the means approach the values in (3.1) and the variances in the runs decrease in the number of historical observations p. Accounting for a potential constraint on the inflow, we compare a scenario where the inflow is limited to 2 (medium constraint) and 1.5 (low constraint). We repeat the idea of Fig.5in Fig.6emphasizing that, except for the constraint, all other quantities remain unchanged. However, the demand illustration seems to be less fluctuating which can be explained by the coarser discretization grid that is used for the constrained optimization. In the medium constraint case, we observe that the inand outflow follow the unconstrained case but are truncated at the very highest peaks and otherwise follow the averaged demand well. Regarding the reconstruction of the weights of the base demand levels when zooming in, one can still observe a quite good match in the inand outflows of the optimized and initial choices of β. Table 2underlines this observation, 123 Inverse demand tracking in transportation networks 655 Fig. 4 Three of the perturbed historically observed demands for demand vertex v7in the case of Sect.3.3.1 Fig. 5 A comparison between the mean realization of the p=6 historical inand outflows with the inand outflow for the reconstructed βin the case of Sect. 3.3.1 and the initial β Table 1 Means and variances of the reobtained weights for the base demands for different choices of the number of perturbed historical observations pin the setting of Sect. 3.3.1 mean variance p=1p=6p=20 p=200 p=1p=6p=20 p=200 β10.2003 0.2000 0.1999 0.2000 2.16e−06 0.25e−06 0.16e−06 0.10e−07 β20.1496 0.1499 0.1500 0.1500 6.65e−06 0.42e−06 0.29e−06 0.21e−07 β30.2001 0.2001 0.2002 0.2000 2.82e−06 0.87e−06 0.31e−06 0.18e−07 β40.4503 0.4500 0.4499 0.4500 4.16e−06 0.57e−06 0.34e−06 0.21e−07 123 656 S. Göttlich et al. Fig. 6 A comparison between the mean realization of the p=6 historical inand outflows with the inand outflow for the reconstructed βin the case of Sect. 3.3.1 and the initial βwith two different inflow constraints Table 2 Means and variances of the reobtained weights for the base demands for different choices of the number of perturbed historical observations pin the setting of Sect.3.3.1 with additional inflow constraints mean variance medium low medium low p=6p=200 p=6p=200 p=6p=200 p=6p=200 β10.1960 0.1972 0.2460 0.2474 2.87e−06 0.11e−06 0.83e−04 0.11e−05 β20.1549 0.1536 0.0997 0.0988 5.93e−06 0.22e−06 0.91e−04 0.10e−05 β30.1997 0.1997 0.1297 0.1284 6.77e−06 0.35e−06 1.10e−04 0.21e−05 β40.4493 0.4495 0.5245 0.5254 6.24e−06 0.30e−06 1.11e−04 0.15e−05 but shows a small deviation especially in the parameters β1and β2compared to the unrestricted case. For the low constraint, the inflow is cut from Monday to Friday and in some peak times also during the weekend, so that most of the time demand cannot be satisfied on average. Then the reconstruction task is also not successful, and we can observe a visible mismatch in the green circles (associated to the optimal outflow for the initial β) and purple diamonds (representing the outflow for the reconstructed β) during the weekend. Referring again to Table 2, one can see that there is a large deviation in the reconstructed values of β, where the very low values of β2and β3 are particularly striking. This effect can be explained by the fact that D2and D3are the sinusoidal components of demand and that the observations are smoothed and truncated at the majority of time. 123 Inverse demand tracking in transportation networks 657 Fig. 7 Three of the perturbed historically observed demands for demand vertex v7in the case of Sect.3.3.2 3.3.2 Results with additional noise in the weights ˇ In addition to the investigation of Sect.3.3.1, we introduce a structural and uncertain deviation in the choice of β, when generating the historically desired demand in (3.2). We assume that the uncertainty mainly comes into play for β4such that for any historical observation, the weights for the demand levels are chosen as (β1,β 2,β 3,β 4)=0.2 1+˜ Z,0.15 1+˜ Z,0.2 1+˜ Z,0.45 +˜ Z 1+˜ Z(3.3) for a uniformly distributed random variable ˜ Z∼U([−0.05,0.05]). The results for some historical observations are presented below in Fig.7. There is not only noise in the demands but also structurally different behavior due to different realizations of ˜ Zin the weights of the demands. Therefore, the yellow curve of historic data 3 seems to be lower (corresponding to a larger value of ˜ Z) than the blue curve (corresponding to a smaller value of ˜ Z). Figure8shows the different inand outflows which are supplemented by Table 3showing the means and the variances of a Monte Carlo simulation for the reconstructed weights of the base demands for different numbers of perturbed historical observations. We observe that in Fig.8,the expected outflow and inflow match quite well, but considering Table 3, it can be seen that the reconstruction is more difficult than in the standard setting. For small p,the reconstructed βdeviates more significantly from the initial choice. For a larger number of observations p, the data indicates that the performances are improved and lead to good reconstructed values of β. Also in this scenario, we investigate a constraint on the inflow control on a medium level of 2 and a low constraint of 1.5. Similar to Sect.3.3.1, the reconstruction works at 123 658 S. Göttlich et al. Fig. 8 A comparison between the mean realization of the p=6 historical inand outflows with the inand outflow for the reconstructed βin the case of Sect. 3.3.2 and the initial β Table 3 Means and variances of the reobtained weights for the base demands for different choices of the number of perturbed historical observations pin the setting of Sect. 3.3.2 mean variance p=1p=6p=20 p=200 p=1p=6p=20 p=200 β10.2012 0.2010 0.2001 0.2001 0.36e−04 0.08e−04 0.19e−05 0.01e−05 β20.1497 0.1509 0.1500 0.1501 0.26e−04 0.04e−04 0.14e−05 0.01e−05 β30.2013 0.2013 0.2003 0.2001 0.40e−04 0.06e−04 0.19e−05 0.01e−05 β40.4477 0.4468 0.4496 0.4497 2.76e−04 0.49e−04 1.37e−05 0.08e−05 least satisfactorily in the medium constraint case, whereas it fails in the low constraint case. Nevertheless, in both cases, the average outflow matches the optimal outflow for the reconstructed β, see Fig.9. Table 4shows for p∈{6,200}the mean and the variance as the adapted version of Table 3with medium and low inflow constraint, where the variances are similar but slightly higher than in the unconstrained framework. The observed effects are comparable to those obtained for the constrained but unperturbed regime in Table 2. 3.3.3 Results with changed base demand level D4 This section is based on the investigations in Sect.3.3.1. Instead of perturbing β,we assume that there is a structural deviation in the base demand levels. Particularly, we assume that in the generation of the observations, we adjust the base demand D4to D4(t)=3 21[0,120](t), which means that there is larger share of demand on weekdays. Furthermore, we omit the normalization restriction to the weights, i.e., we merely assume β≥0, ∈{1,...,4}, and drop the constraint 4 =1β=1, since the increase in the base demand level should now be captured by a larger weight on 123 Inverse demand tracking in transportation networks 665 affine constraints, see the classical paper (Dempe and Bard 1992) for a related idea and Gfrerer and Outrata (2024) for a modern view. A Special quadratic problems with box constraints Let us fix vectors θ,vd∈Rnas well as va∈(R∪{−∞})nand vb∈(R∪{∞})nsuch that all entries of θare positive while va≤vbholds componentwise. For := diag(θ), we aim to solve min v{1 2vv −v dv|v∈Vad}(QP) where Vad ⊂Rnis the box given by Vad := {v∈Rn|va≤v≤vb}. First, we observe that the objective function in (QP) is uniformly convex while the feasible set is nonempty, closed, and convex. Hence, (QP) possesses a uniquely determined global minimizer ¯v∈Vad. The latter can be characterized in terms of the necessary and sufficient optimality condition ∀v∈Vad :(¯v−vd)(v −¯v) ≥0.(A.1) We note that is a positive definite diagonal matrix. Hence, it is reasonable to set ˜v:= −1vd. Note that ˜vi=θ−1 ivd,iholds for all i=1,...,n. We will now show that ¯v=max(va,min(˜v,vb)) (A.2) holds true, i.e., that ¯vis the projection of ˜vonto the box Vad. Note that max and min have to be interpreted componentwise in (A.2). We introduce index sets Ia,I0,Ib⊂ {1,...,n}by means of Ia:= {i∈{1,...,n}|˜vi<v a,i}, I0:= {i∈{1,...,n}|va,i≤˜vi≤vb,i}, Ia:= {i∈{1,...,n}|vb,i<˜vi}. Clearly, these sets form a disjoint partition of {1,...,n}, and (A.2) can be rewritten as ∀i∈{1,...,n}: ¯vi=⎧ ⎪ ⎨ ⎪ ⎩ va,ii∈Ia, ˜vii∈I0, vb,ii∈Ib. 123 666 S. Göttlich et al. Pick v∈Vad arbitrarily. Taking together all of the above findings, we end up with (¯v−vd)(v −¯v) = i∈Ia (θiva,i−vd,i)(v i−va,i)   ≥0 + i∈Ib (θivb,i−vd,i)(v i−vb,i)   ≤0 + i∈I0 (θi˜vi−vd,i)   =0 (vi−˜vi) ≥ i∈Ia (θi˜vi−vd,i)   =0 (vi−va,i)+ i∈Ib (θi˜vi−vd,i)   =0 (vi−vb,i)=0, and this shows that ¯vconstructed as in (A.2) is, indeed, a solution of (A.1) and, thus, the uniquely determined global minimizer of (QP). Acknowledgements The authors wish to thank the two anonymous reviewers whose valuable comments and suggestions helped to improve the overall quality of this paper. Furthermore, one of the reviewers recommended an inspection of the PhD thesis (Keimer 2014) which is gratefully acknowledged. Simone Göttlich was supported by the Deutsche Forschungsgemeinschaft (DFG) within the projects GO1920/10-1 and GO1920/11-1. Funding Open Access funding enabled and organized by Projekt DEAL. Declarations Conflict of interest The authors declare no Conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References Albrecht S, Ulbrich M (2017) Mathematical programs with complementarity constraints in the context of inverse optimal control for locomotion. Optim Methods Softw 32(4):670–698. https://doi.org/10. 1080/10556788.2016.1225212 Albrecht S, Passenberg C, Sobotka M, Peer A, Buss M, Ulbrich M (2010) Optimization criteria for human trajectory formation in dynamic virtual environments. In: Kappers AML, van Erp JBL, Bergmann Tiest WM, van der Helm FDT (eds) Haptics: generating and perceiving tangible sensations. Springer, Berlin, pp 257–262. https://doi.org/10.1007/978-3-642-14075-4_37 Albrecht S, Leibold M, Ulbrich M (2012) A bilevel optimization approach to obtain optimal cost functions for human arm movements. Numer Algebra Control Optim 2(1):105–127. https://doi.org/10.3934/ naco.2012.2.105 Banda MK, Herty M, Klar A (2006) Gas flow in pipeline networks. Netw Heterog Media 1(1):41–56. https://doi.org/10.3934/nhm.2006.1.41 Bard JF (1998) Practical bilevel optimization. Springer, New York. https://doi.org/10.1007/978-1-47572836-1 Bressan A (2000) Hyperbolic systems of conservation laws-the one-dimensional cauchy problem. Oxford University Press, Oxford. https://doi.org/10.1093/oso/9780198507000.001.0001 Bressan A, ˇ Cani´c S, Garavello M, Herty M, Piccoli B (2014) Flows on networks: recent results and perspectives. EMS Surv Math Sci 1(1):47–111. https://doi.org/10.4171/EMSS/2 123 Inverse demand tracking in transportation networks 667 Clarke FH (1983) Optimization and nonsmooth analysis. Wiley, New York. https://doi.org/10.1137/1. 9781611971309 Colton D, Kress R (2013) Inverse acoustic and electromagnetic scattering theory. Springer, New York. https://doi.org/10.1007/978-1-4614-4942-3 Coskun S, Korn R (2021) Modeling the intraday electricity demand in Germany. In: Göttlich S, Herty M, Milde A (eds) Mathematical modeling, simulation and optimization for power engineering and management. Springer, Cham, pp 3–23. https://doi.org/10.1007/978-3-030-62732-4_1 Dempe S (2002) Foundations of bilevel programming. Kluwer, Dordrecht. https://doi.org/10.1007/b101970 Dempe S (2020) Bilevel optimization: theory, algorithms, applications and a bibliography. In: Dempe S, Zemkoho AB (eds) Bilevel optimization: advances and next challenges. Springer, Cham, pp 581–672. https://doi.org/10.1007/978-3-030-52119-6_20 Dempe S, Bard JF (1992) Bundle trust-region algorithm for bilinear bilevel programming. J Optim Theory Appl 110:265–288. https://doi.org/10.1023/A:1017571111854 Dempe S, Kalashnikov V, Pérez-Valdéz G, Kalashnykova N (2015) Bilevel programming problems-theory, algorithms and applications to energy networks. Springer, Berlin. https://doi.org/10.1007/978-3-66245827-3 Dempe S, Harder F, Mehlitz P, Wachsmuth G (2019) Solving inverse optimal control problems via value functions to global optimality. J Global Optim 74(2):297–325. https://doi.org/10.1007/s10898-01900758-1 Dobrowolski M (2006) Angewandte Funktionalanalysis. Springer, Berlin. https://doi.org/10.1007/3-54029960-2 Friedemann M, Harder F, Wachsmuth G (2023) Finding global solutions of some inverse optimal control problems using penalization and semismooth Newton methods. J Global Optim 86:1025–1061. https:// doi.org/10.1007/s10898-023-01288-7 Gfrerer H, Outrata JV (2024) On the role of of semismoothness in nonsmooth numerical analysis: theory. https://arxiv.org/abs/2405.14637 Göttlich S, Schillinger T (2022a) Control strategies for transport networks under demand uncertainty. Adv Comput Math 48:74. https://doi.org/10.1007/s10444-022-09993-9 Göttlich, S, Schillinger T (2022b) Stochastic optimal control for nonlinear damped network dynamics. arXiv:2202.05114 Göttlich S, Herty M, Schillen P (2016) Electric transmission lines: control and numerical discretization. Optimal Control Appl Methods 37(5):980–995. https://doi.org/10.1002/oca.2219 Göttlich S, Korn R, Lux K (2019) Optimal control of electricity input given an uncertain demand. Math Methods Oper Res 90:301–328. https://doi.org/10.1007/s00186-019-00678-6 Gugat M, Keimer A, Leugering G, Wang Z (2015) Analysis of a system of nonlocal conservation laws for multi-commodity flow on networks. Netw Heterog Media 10(4):749–785. https://doi.org/10.3934/ nhm.2015.10.749 Gugat M, Schultz R, Wintergerst D (2018) Networks of pipelines for gas with nonconstant compressibility factor: stationary states. Comput Appl Math 37(2):1066–1097. https://doi.org/10.1007/s40314-0160383-z Harder F, Wachsmuth G (2019) Optimality conditions for a class of inverse optimal control problems with partial differential equations. Optimization 68(2–3):615–643. https://doi.org/10.1080/02331934. 2018.1495205 Hatz K, Schlöder JP, Bock HG (2012) Estimating parameters in optimal control problems. SIAM J Sci Comput 34(3):A1707–A1728. https://doi.org/10.1137/110823390 Hinze M, Pinnau R, Ulbrich M, Ulbrich S (2009) Optimization with PDE constraints. Springer, Dordrecht. https://doi.org/10.1007/978-1-4020-8839-1 Holler G, Kunisch K, Barnard RC (2018) A bilevel approach for parameter learning in inverse problems. Inverse Prob 34(11):1–28. https://doi.org/10.1088/1361-6420/aade77 Keimer A (2014) Optimal control of nonlinear nonlocal conservation laws on networks. PhD thesis, University of Erlangen–Nuremberg. https://open.fau.de/items/b9a127b7-0e85-43d0-822b-f391b6f206ea Mehlitz P, Wachsmuth G (2020) Bilevel optimal control: existence results and stationarity conditions. In: Dempe S, Zemkoho AB (eds) Bilevel optimization: advances and next challenges. Springer, Cham, pp 451–484. https://doi.org/10.1007/978-3-030-52119-6_16 Mombaur K, Truong A, Laumond J-P (2010) From human to humanoid locomotion–an inverse optimal control approach. Auton Robot 28(3):369–383. https://doi.org/10.1007/s10514-009-9170-7 123 668 S. Göttlich et al. Rein M, Mohring J, Damm T, Klar A (2020) Optimal control of district heating networks using a reduced order model. Optimal Control Appl Methods 41(4):1352–1370. https://doi.org/10.1002/oca.2610 Schramm H, Zowe J (1992) A version of the bundle idea for minimizing a nonsmooth function: conceptual idea, convergence analysis, numerical results. SIAM J Optim 2(1):121–152. https://doi.org/10.1137/ 0802008 Shimizu K, Ishizuka Y, Bard JF (1997) Nondifferentiable and two-level mathematical programming. Springer, New York. https://doi.org/10.1007/978-1-4615-6305-1 Sikolya E (2004) Semigroups for flows in networks. PhD thesis, University of Tübingen. https://d-nb.info/ 972906657/34 Suryan V, Sinha A, Malo P, Deb K (2016) Handling inverse optimal control problems using evolutionary bilevel optimization. In: 2016 IEEE congress on evolutionary computation (CEC), pp 1893–1900. https://doi.org/10.1109/CEC.2016.7744019 Tröltzsch F (2010) Optimal control of partial differential equations. Am Math Soc. https://doi.org/10.1090/ gsm/112 Troutman JL (1996) Variational calculus and optimal control. Springer, New York. https://doi.org/10.1007/ 978-1-4612-0737-5 Vinter R (2010) Optimal control. Birkhäuser, Boston. https://doi.org/10.1007/978-0-8176-8086-2 Zemkoho AB (2016) Solving ill-posed bilevel programs. Set-Valued Var Anal 24:423–448. https://doi.org/ 10.1007/s11228-016-0371-x Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123