FREE BOUNDARY VALUE PROBLEMS AND HJB EQUATIONS FOR THE STOCHASTIC OPTIMAL CONTROL OF ELASTO-PLASTIC OSCILLATORS

Abstract. We consider the optimal stopping and optimal control problems related to stochastic variational inequalities modeling elasto-plastic oscillators subject to random forcing. We formally derive the corresponding free boundary value problems and Hamilton-Jacobi-Bellman equations which belong to a class of nonlinear partial of differential equations with nonlocal Dirichlet boundary conditions. Then, we focus on solving numerically these equations by employing a combination of Howard’s algorithm and the numerical approach [A backward Kolmogorov equation approach to compute means, moments and correlations of non-smooth stochastic dynamical systems; Mertz, Stadler, Wylie; 2017] for this type of boundary conditions. Numerical experiments are given.


Motivations and background
In an extremely wide range of applications, mechanical systems are fundamentally affected by vibrations. For instance, components (such as pipes) in a power plant are designed to be structurally robust under seismic vibrations. Vibrations also cause mechanical systems to accumulate plastic deformations and eventually fail. The importance of these issues has motivated a lot of effort in modeling. From the modeling point of view, keeping track of the impact of past vibrations requires a specific description of the mechanical systems under consideration. The state of the system must be described by a randomly forced dynamical system with memory. Random forcing here expresses the stochastic nature of the vibrations that apply to the mechanical structures. In this framework, the difficulty is to handle dynamical systems with memory. A huge engineering literature has been devoted to this topic [14, 17, 18, 25, 29, 31, 34-36, 40, 41]. This paper mainly focuses on an important model, referred to as the elastic-perfectly-plastic oscillator (EPPO), which was introduced in [28] and appears in various applications [16, 19-21, 23, 24]. In terms of rheology, the response of such model is a combinations of elementary rheological models : springs, Saint-Venant's elements, dashpots and material points as shown in Figure 1.
η Figure 1. Rheological model of an elasto-perfectly-plastic oscillator. A mass (black box) is associated in series with elements which are themselves an association in parallel or in series of elementary rheological models. The whole system is excited by a time-dependent force η.
In [12,13], the stochastic variational inequalities (SVIs) have been identified as the right mathematical tool to deal with this model. Here, our objective is to derive and solve numerically a new class of free boundary problems and Hamilton-Jacobi-Bellman equations related to the control of these SVIs. Below, we provide background on and identify research issues on a stochastically driven EPPO.

An elasto-perfectly-plastic problem with noise
The dynamics focuses on two quantities: the deformation supported by the structure, namely x t , when subjected to the vibrations, and its velocity, namely y t " 9 x t . In the simplest case, the dynamics of the velocity is given by 9 y t`c0 y t`ft " η t , where c 0 ą 0 is a damping coefficient, η t describes the random vibrations and f t describes the restoring force arising from the structure. The exact nature of f t depends on the particular structure under consideration. The so-called perfect model corresponds to the case in which the force f t is a linear function of x t in the elastic phase and a constant in the plastic phase. In this case, the permanent (plastic) deformation at time t can be written as ∆ t fi ş t 0 1 t|fs|"kY u y s ds, where kY describes the plastic phase through a stiffness coefficient k and an elasto-plastic bound, that is normalized to Y " 1. The variable z t " x t´∆t describes the deformation without plastic deformation. When |z t | ă 1, the structure is in the elastic regime and when |z t | " 1, it is in the plastic regime. In both cases, the restoring force f t is given by kz t . The dynamics is then described by the pair py t , z t q that satisfies 9 y t`c0 y t`k z t " η t , p 9 z t´yt qpφ´z t q ě 0, @|φ| ď 1, |z t | ď 1, @t.

Research Issues
From a mathematical perspective, the point is to understand the behavior of the model in order to predict the behavior of the mechanical structure of interest. For engineering applications, this is fundamental in assessing the risk of failure of the structure. In [5,6,12], the ergodic behavior of the couple velocity-restoring force py t , z t q of an elasto-perfectly-plastic oscillator (EPPO) excited by white noise and by a filtered noise has been established. The invariant probability can be described, through its density, as the solution of a partial differential equation (PDE), with specific boundary conditions arising from the plastic properties of the model. It is worth mentioning that ergodicity does not hold for the plastic deformation ∆ t . Thus, the behavior of ∆ t can be investigated through different quantities such as (growth rate) variance of the plastic deformation E " ∆ 2 T ‰ T and the failure risk Pˆmax where b is a critical threshold of deformation considered as dangerous for the system and p0, T q is a given interval of time. For large time, analysis of the variance for the plastic deformation (2) has been carried out in [8] for an EPPO excited by white noise and a probabilistic formula has been derived. Also, in [7] an analytical formula has been proven. Concerning the failure risk, the mathematical techniques that have been developed so far provide an asymptotic formula [30]. Analysis of the risk raises other important problems for which no answer has been given so far. ‚ The risk probability (2) can be reformulated in terms of an optimal stopping problem as follows: This reformulation relies on a general result on reachability for Markov processes [10]. The problem is technically delicate 1 . Here, we rather focus on a similar problem for the plastic state, that is, In subsection 2.1, we describe a general framework based on a new class of free boundary value problems characterizing optimal stopping problems of the form sup τ PT0,T Eˆż τ 0 e´λ s gpz s , y s qds`e´λ τ f pz τ , y τ q˙.
‚ Another important concept related to the notion of critical excitation [38] must be considered. The major point is to answer the following question: can we find the most likely random force, of the form η t " α t`9 w t where w t is a Wiener process, that would maximize the accumulation of the plastic deformation or the time spent in plastic regime or in the sense sup αp.q We aim at translating the problem within the framework of stochastic control. In subsection 2.2, we describe a general framework based on a new class of Hamilton-Jacobi-Bellman Equations aimed at characterizing the control α in the form sup αp.q E˜ż T 0 e´λ s gpz s , y s qds`f pz T , y T q¸.

Notations
Given T ą 0, set D T fi p´1, 1qˆRˆr0, T s, DT fi t1uˆRˆr0, T s and DT fi t´1uˆRˆr0, T s. Then, define three differential operators A, B`and B´as follows Bz .
We use the notation C ‹ T for the set of continuous functions on r´1, 1sˆRˆr0, T s that are C 1 -regular with respect to z, C 2 -regular with respect to y and C 1 -regular with respect to t.

Overview of our numerical study
The PDEs related to the optimal stopping and optimal stochastic control of the elasto-plastic oscillator are nonstandard free boundary value and HJB problems. To the best of our knowledge, it seems they have not been treated in the literature. In the linear case that is without control (i.e backward Kolmogorov equations), only partial existence and uniqueness results are available [9,13]. This is because standard PDE theory techniques do not apply due to the non-standard boundary conditions and the degeneracy of these problems. Therefore, in the note, we mostly study the behavior of these PDEs numerically in order to gather insight in the solution behavior.

Derivation of the free boundary value and HJB problems
In this section, we derive formally a new class of free boundary value problems for the optimal stopping (see subsection 2.1) and a new class of Hamilton-Jacobi-Bellman Equations for the optimal control (see subsection 2.2) of the elasto-plastic oscillator (1). These two problems are time dependent on a finite time interval with a terminal condition. Infinite horizon and stationary versions of these two problems are also given (see subsection 2.3).
2.1. Free boundary value problem for the optimal stopping of (1) We are not aware of mathematical studies regarding optimal stopping problems for SVIs modeling elasto-plastic oscillators and their applications to the risk analysis of failure. Thus, to address this issue, the strategy here is to rely on the connection between optimal stopping and variational inequalities as established by Bensoussan-Lions [3]. For an EPPO excited by white noise, this approach leads to a new type of free boundary value problem governed by a variational inequality with a nonlocal Dirichlet boundary condition.
Theorem 1 (a free boundary value problem with nonlocal Dirichlet boundary condition for an EPPO excited by a white noise). Let λ P R and u P C ‹ T , assume that u satisfies maxˆB u Bt`A u´λu`g, f´u˙" 0 in D T and maxˆB u Bt`B˘u´λ u`g, f´u˙" 0 in DT with the terminal condition upz, y; T q " f pz, yq. Then upz, y, tq has a probabilistic interpretation in terms of an optimal stopping problem, that is upz, y, tq " sup τ PTt,T Eˆż τ t e´λ ps´tq gpz s , y s qds`e´λ pτ´tq f pz τ , y τ q|pz t , y t q " pz, yqẇ here pz s , y s q satisfies (1) with η is a white noise and T t,T is the set of stopping times between t and T .
Proof. Without loss of generality, we can assume t " 0 and prove upz, y, 0q " sup where pz s , y s q satisfies (1). Define We claim that M t is a martingale. Since τ P T 0,T is bounded, by optional sampling theorem, E rM τ s " E rM 0 s . Also using the inequalities in the hypothesis on u, we have for any t P r0, T s, M t ě e´λ t upz t , y t , tqş t 0 e´λ s gpz s , y s qds. So on one hand, upz, y, 0q From the terminal condition of u, we know τ opt P T 0,T . Meanwhile for any s P r0, τ opt q,upz s , y s , sq ą f pz s , y s q, so from the hypothesis we must have Bu Bt`A u´λu`g " 0 on D T and Bu Bt`B˘u´λ u`g " 0 on DT . Therefore using optional sampling theorem again and the continuity of u, we have Collecting results, the probabilistic interpretation is obtained.

HJB equations for the stochastic control problem
A natural approach to address the problem of the critical excitation [38] is the Bellman principle.

Statement of the problem
Let 0 ď t ď T , α : rt, T s Ñ R be a given function and consider the state variable pz α , y α q satisfying (1) with η " α`9 w, that is 9 y α t`c0 y α t`k z α Remark 1. By changing the probability measure, which is indeed possible by means of Girsanov transform, we can easily get rid of the control α (say when it is bounded). This permits to address the solvability of the SVI in the weak sense, provided that the underlying filtration is chosen accordingly. This technique is standard for SDEs and for the so-called "weak formulation" of related optimal stochastic control problems.
In this context, we consider α as a control variable for the state variable pz α , y α q. Next, we define the functional where pz α s , y α s q sět satisfies (8) with pz α s , y α s q " pz, yq. For instance, if λ " 0, f " 0 and gpz, yq " |y|1 |z|"1 then the integrand in J represents the accumulated plastic deformation. We say that the noiseα z,y,t`9 w is a critical excitation whenα z,y,t is an optimal control for Problem (8)-(9) maximizing the functional J z,y,t , that is, sup αp¨qPAt,T J z,y,t pαp.qq " J z,y,t pαp.qq. Here A t,T is a set of functions taking values in A " r´1, 1s. The problem is to find the value function vpz, y, tq " sup αp¨qPAt,T J z,y,t pαp.qq and an optimal controlα z,y,t .

Bellman principle and HJB equation
Bellman principle can be formulated in our context as follows: if we consider an optimal controlα |rt,T s for Problem (8)-(9) starting at time t with initial condition pz, yq, then for any intermediate time t`h between the times t and T , the portion of control enclosed by the times t`h and T ,α |rt`h,T s remains optimal for the problem starting at time t`h with initial condition pzα t`h , yα t`h q resulting fromα |rt,t`hs . Mathematically, it reads The infinitesimal version of the Bellman principle leads to a new type of HJB equation with a nonlocal Dirichlet boundary condition.
Theorem 2 (A nonlocal HJB problem related to an EPPO).
Let v P C ‹ T and assume that it satisfies with the terminal condition vpz, y, T q " f pz, yq and where, for any α P A, A α fi α B By`A , B α fi α B By`B˘.
Moreover, we assume that (a) there existsâpz, y, tq maximizing a Þ Ñ`B v Bt`A a v´λv`g˘pz, y, tq andâp˘1, y, tq tq is a well-defined control process in A and (c) the SVI in (8) with α t "α t defines a unique solution pẑ s ,ŷ s q sět for each given initial data pẑ t ,ŷ t q " pz, yq. Then vpz, y, tq has the following probabilistic interpretation as the value function of an optimal stochastic control problem, vpz, y, tq " sup αPA J z,y,t pαp.qq.
Going back to the case where λ " 0, f " 0 and gpz, yq " |y|1 |z|"1 , the objective functional so findingα would allow us to know a certain type of random forcing that would maximize the accumulated plastic deformation over the interval of time p0, T q.
Proof. In our context, the proof follows the line of [39]. Let αp.q P A t,T be an arbitrary control, pz α s , y α s q sět the associated state process with initial condition pz t , y t q " pz, yq. By Itô's lemma applied to v between t and T , But since for every s P rt, T s, Bv Finally, from the (a)-(b)-(c) hypotheses, the control achieves the equality in (12).

Infinite horizon problems
A similar set of problems as those presented in equations (5) and (10) can be obtained in the infinite time horizon and stationary case. We give the two statements without proof since they are similar to what was done for the time dependent case. As the time horizon is infinite, there are a few additional assumptions.
Theorem 3 (an infinite horizon free boundary value problem with nonlocal Dirichlet boundary condition for an EPPO excited by a white noise). Let λ ě 0 and u P C ‹ T "8 , assume that u satisfies max pAu´λu`g, f´uq " 0 in D and max pB˘u´λu`g, f´uq " 0 in D˘.
Then upz, yq has a probabilistic interpretation in terms of an optimal stopping problem, that is upz, yq " sup τ PT Eˆż τ 0 e´λ s gpz s , y s qds`e´λ τ f pz τ , y τ q|pz t , y t q " pz, yqẇ here pz s , y s q satisfies (1) with η is a white noise and T is the set of stopping times almost surely finite.
Theorem 4 (an infinite horizon nonlocal HJB problem related to an EPPO). Let λ ě 0 and v P C ‹ T "8 , assume that it satisfies max αPA tA α v´λv`gu " 0 in D and max αPA B α v´λv`g ( " 0 in D˘ (14) where, for any α P A, A α fi α B By`A , B α fi α B By`B˘. Moreover, we assume that paq 8 there existsâpz, yq maximizing a Þ Ñ pA a v´λv`gq pz, yq andâp˘1, yq maximizing a Þ Ñ`B ȃ v´λv`g˘p˘1, yq such that Aâ pz,yq v´λv`g¯pz, yq " 0 in D, and´Bâ p˘1,yq v´λv`g¯p˘1, yq " 0 in D˘, pbq 8 the processα t fiâpz t , y t q is a well-defined control process in A and pcq 8 the SVI in (8) with α t "α t defines a unique solution pẑ s ,ŷ s q sět for each given initial data pẑ t ,ŷ t q " pz, yq. Then vpz, yq has the following probabilistic interpretation as the value function of an infinite horizon optimal stochastic control problem, vpz, yq " sup αPA J z,y,t pαp.qq where J z,y pαp¨qq fi E "ż 8 0 e´λ s gpz α s , y α s qdsˇˇpz 0 , y 0 q " pz, yq Let us comment on the shape of the optimizerâ. Since g does not depend on α, Equation (14) can be simply written as max αPA " α Bv By Since A " r´1, 1s, we must have nd the maximizer is given byâpz, yq " sign´B v By pz, yq¯, where signpxq fi x |x| , if x ‰ 0 and 0 if x " 0. It turns out that the (feedback) controlα t "âpz t , y t q takes values in t´1, 1u and thus belongs to the class of bang-bang controls.
Furthermore, when f, g are symmetric functions i.e f pz, yq " f p´z,´yq and gpz, yq " gp´z,´yq, it can be seen thatâ must be an antisymmetric function with respect to py, zq. Indeed, denotingṽ : pz, yq Þ Ñ vp´z,´yq, it can be readily checked thatṽ solves the same PDE as v (using the fact that A is a symmetric interval around 0). Assuming uniqueness of the solution to this PDE, we haveṽ " v. Hence v is symmetric and thus Bv By andâ are antisymmetric. A similar discussion holds for the time dependent problem in Theorem 2.
Numerically, in the case of accumulated plastic deformation (see Figures 8-9 below), we have indeed found an antisymmetric function that seems to take the value 1 (resp.´1) on a set of the form tpx, yq, γpxq ą yu (resp. tpx, yq, γpxq ă yu). The interface γ seems to be of the form γpxq " cx where c is a certain constant. The study of this interface is an interesting question that we leave for future work.

Howard's algorithm
In this section, we recall Howard's algorithm [1,2,26] which treats the following typical problem: find v P R N , satisfying the componentwise minimization where A is a nonempty compact set and for all α P A N , Bpαq is a NˆN matrix and cpαq is a vector of size N . As a particular case, it also includes the following problem: find v P R N , satisfying the componentwise minimization min´Bv´g , v´f¯" 0 whereB is a NˆN matrix and h, g are vectors of size N . Indeed, it can be written under the form (16) with, for all 1 ď i ď N , HereB i,‚ is the line i of the matrixB and I is the NˆN identity matrix. As expressed for example in [11], Howard's algorithm for (16) is as follows: (1) Initialize α p0q P A N (2) Iterate for k ě 0: (i) Find v pkq P R N solution of Bpα pkq qv pkq´c pα pkq q " 0. If k ě 1 and v pkq " v pk´1q then stop. Otherwise go to (ii).
(ii) α pk`1q fi arg min aPA N`B paqv pkq´c paq˘. Set k fi k`1 and go to (i).
(ii) For every 1 ď i ď N, α pk`1q i fi # 0 if pBv pkq´g q i ď pv pkq´f q i 1 otherwise. Set k fi k`1 and go to (i). In the numerical implementation, we replace the stopping condition v pkq " v pk´1q by a criterion of the form For convergence results related to this algorithm, we refer the reader to [11]. Now we explain how to solve numerically equations (5), (10), (13) and (14) using this algorithm. To this end, we need to discretize these equations and rewrite them under the form (16) and (17).

Discretization of PDEs problems following [32]
To numerically approximate the solutions of (5), (10), (13) and (14), we use a finite difference scheme. We truncate the unbounded domains D and D˘to obtain D Y fi p´1, 1qˆp´Y, Y q and DY fi t˘1uˆp´Y, Y q, where Y is chosen sufficiently large that the probability of finding the underlying process outside D Y and DY is negligible. We apply a homogeneous Neumann boundary condition at y "˘Y . We consider a two-dimensional rectangular finite difference grid, G fi tpz i , y j q fi p´1`pi´1qδz,´Y`pj´1qδyqu 1ďiďI,1ďjďJ , where δz fi 2 I´1 , δy fi 2Y J´1 . Here, I, J are odd integers of the form 2Ĩ`1, 2J`1. The total number of nodes in G is N " IJ. The numerical approximations of upz i , y j , t n q, vpz i , y j , t n q, u λ pz i , y j q and v λ pz i , y j q are denoted by u n i,j , v n i,j , u λ i,j and v λ i,j and the corresponding vectors collecting all the unknowns are u n , v n , u λ and v λ . We use the notation f i,j , g i,j for f px i , y j q, gpx i , y j q and the corresponding vectors are f , g. Here, t n fi nδt discretizes the time and N T δt " T . (5) and (13) For the time dependent problem (5), we use an implicit Euler method to discretize in time together with finite differences in space and, for the stationary problem (13), we only use finite differences in space. For both problems, the condition on D in (5) and in (13) at the black points in Figure 2 result in u 0 i,j " f i,j and min˜u n`1 i,j´u n i,j δt´`L u n`1˘i ,j´g i,j , u n`1 i,j´f i,j¸" 0, 0 ď n ď N T´1 ,

Case of the free boundary problems
and min´λu λ i,j´`L u λ˘i ,j´g i,j , u λ i,j´fi,j¯" 0, where pLuq i,j is of the forḿ pLuq i,j fi´pL y uq i,j´m axp0, y j qˆu i`1,j´ui,j δz˙´m inp0, y j qˆu Figure 2. Discretization of D Y " p´1, 1qˆp´Y, Y q and DY " t˘1uˆp´Y, Y q. At black points, the discretized equation is satisfied. At grey points, homogeneous Neumann boundary conditions are used, and at red points non-standard boundary conditions are employed. with´p δyẇ here b i,j "´pc 0 y j`k z i q. The nonlocal Dirichlet boundary conditions (second condition) in (5) and in (13) at the red points in Figure 2 are discretized by the same formulae with L replaced by L˘, where pL˘uq i,j are defined by´p L`uq i,j fi´pL y uq i,j´m inp0, y j qˆu i,j´ui´1,j δz˙ ( 21) and´p L´uq i,j fi´pL y uq i,j´m axp0, y j qˆu i`1,j´ui,j δz˙.
Note that, unlike [11,27], we are not able to guess the value of u and u λ or its derivative on the boundaries y "˘Y . Instead, we proceed as in [32] in imposing an homogeneous Neumann condition. This is based on the intuition that, when Y is chosen sufficiently large, the probability of finding the underlying process outside the truncated domain is negligible. This idea is supported by Monte Carlo simulations as shown below.
The Neumann boundary conditions at the points shown in grey in Figure 2 results in We also do the following modification : g i,j " 0 at the grey points. For the time dependent problem, this results in the following nonlinear system to be solved in each time step: min`M δt u n`1´gn δt , u n`1´f˘" 0, u 0 " f , and for the stationary problem min`M λ u λ´g , u λ´f˘" 0.
The construction of M λ is very similar. For the computational results presented in the remainder of this paper, we use a C code that implements a Monte Carlo (MC) probabilistic simulation to approximate the solution of (1). We use a MATLAB implementation for the PDE approach (24). Implementations are available upon request. (10) and (14) We proceed as what was done above, except that we need to incorporate the control in the upwinding differences. The first conditions in (10) and (14) at the black points in Figure 2
Again, the nonlocal Dirichlet boundary conditions (second condition) in (10) and (14) at the red points in Figure 2 are discretized by the same formulae with L α replaced by L α defined bý pL ὰ uq i,j fi´pL α y uq i,j´m inp0, y j qˆu i,j´ui´1,j δzȧ nd´p L ά uq i,j fi´pL α y uq i,j´m axp0, y j qˆu i`1,j´ui,j δz˙. Also, we impose a Neumann boundary condition at the points shown in grey as in (23). Here, for the time dependent problem, this results in the following nonlinear system to be solved in each time step: and for the stationary problem min where the matrices M α δt and M α λ are almost the same matrices as M δt and M λ , except that b is replaced by b`α.

Free boundary value problems
We consider the case shown in (3) which is related to the problem (5) with λ " 0.1, g " 0 and f " 1 t|z|"1u . This quantity is the probability that the elasto-plastic system has entered, at least once, in the plastic phase before the instant T . Roughly speaking, as long as it remains small, it is somehow reasonable to approximate statistics of the elasto-plastic system by using those of a linear system. By solving (24), a numerical approximation for up0, 0, t n q of (5) is obtained. For comparison, a probabilistic numerical approximation of the left hand side in (3) is provided as an alternative method to the deterministic approach. The numerical results are shown in Figure 3. Note that in contrast with the Monte Carlo approach, the PDE approach also computes upx i , y j , tq, @i, j, as presented in Figure 4 (with λ " 0.1). Figure 5 displays the shape of the optimal stopping decision. We see that, except for the initial time, the optimal stopping decision is constant in time. Recall that, following the notations introduced in Section 3, 1 corresponds to stopping and 0 to continuing. We see that the optimal decision is to stop when the process reaches the plastic boundary and to continue otherwise. Note that the process never reaches the regions where signpyq z "´1 hence the value of the control in these regions does not matter.

HJB problems
We consider both cases shown in (4) which are related to the problem (10) with λ " 0.1, f pz, yq " 0 and ‚ either gpz, yq " |y|1 t|z|"1u (accumulated plastic deformation), ‚ or gpz, yq " 1 t|z|"1u (time spent on the plastic boundary). By solving (28), a numerical approximation for vp0, 0, t n q of (10) is obtained. We do not have explicit solutions for the problems we are studying here and we cannot compute them using a straightforward probabilistic approach. However, for checking, we consider a numerical approximation of an optimal control α opt p.q satisfying  (5) and MC (dashed) approach based on the right hand side of (3). Approximation of the values up0, 0, T q and P 0,0 pDt P r0, T s, |z t | " 1q are shown for T P r0, 7s.
of Monte Carlo simulations based on an Euler discretization of the dynamics (8) controlled by α opt , see e.g. Bernardin [15] for more details. This method has given satisfactory results in the sense that the empirical expectation computed by the probabilistic method matches the corresponding value of u solving (28) computed by the deterministic method. In order to check the convergence of the method, we compute the residual in the sup-norm for (28) with pu opt , α opt q obtained numerically with the above PDE method, that is max i,j,n`p I`δtpM`M αopt qqu n`1 opt´u n opt´δ t g˘i ,j .
For a grid of 400 time steps and I " J " 150, we obtain a residual of the order of 10´1 2 . Figure 6, shows the graph of vp0, 0, T q{T for T P r0, 100s, computed by the PDE method described in Section 3, for the two test cases accumulated plastic deformation (i.e. gpz, yq " |y|1 t|z|"1u ) and time spent on the plastic boundary (i.e. gpz, yq " 1 t|z|"1u ). One can see that for large T , the time spent in plastic deformation is smaller than the accumulated plastic deformation. This is because for large time, even when it starts from p0, 0q, the process has a high probability to reach large values of |y| on the boundary |z| " 1 and when |y| ą 1, the running cost of the test case time spent on the plastic boundary is smaller than the running cost of the other test case. Figure 7 shows the time evolution of the value function for the accumulated plastic deformation test case (with λ " 0.1). The range of the value function keeps increasing with time because with longer time horizon, the process will accumulate more plastic deformation. However, for a given z, the value function is increasing with respect to |y|. This can be explained as follows: when the process starts closer to the plastic boundary, it is more likely to accumulate a lot of plastic deformation, which is reflected by the value function of this problem. The shape optimal control (at the final time t " 100) is displayed in Figure 8 (3d plot) and Figure 9 (2d plot). One can see that, roughly speaking, the optimal decision is to push the process towards the closest portion of the plastic boundary. One can see that the control is´1 in a region and`1 in another region. The line between them crosses the quadrant tpy, zq : y ą 0, z ă 0u and tpy, zq : y ă 0, z ą 0u and, as t increases, it becomes closer to the line y " 0. The exact evolution of this line is an interesting question that we leave for future work. Using the optimal control computed from this PDE method, we can then simulate trajectories of the optimally controlled process. Figure 10 shows, for one realization of the noise, the trajectory with the optimal control (for the sake of completeness, we also show the non-controlled trajectory), and the   (5) : evolution of the value function in time upy, z, tq. The shots are taken at different times t " 0, 1, 5, 10. The terminal condition is f " 1 t|z|"1u , the right hand side g " 0 and λ " 0.1.
corresponding evolution of the optimal control. One can see that the value of the optimal control switches between´1 and`1. Moreover, the value`1 is applied, roughly speaking, in the upper half of the phase space, while the value´1 is applied in the lower part. This is consistent with Figures 8-9. The value function and the optimal control for the problem of the time spent in the plastic deformation are very similar to those obtained for the case of accumulated plastic deformation so we omit them. Furthermore, we have also solved the ergodic problems (free boundary value problem and Hamilton-Jacobi-Bellman problem) and checked that the value functions of dynamic problems converge to the one of the corresponding ergodic problems.

Conclusion
In this work we have investigated control problems for elasto-plastic oscillators. We have derived the PDEs for free boundary value problems and Hamilton-Jacobi-Bellman problems. We then implemented numerical methods based on Howard algorithm to compute approximate solutions. From here, several directions of research are still open. First, it would be interesting to clarify the structure of the optimal control for the problem of  (5): evolution of the optimal stopping decisionâpz, y, tq as shown in (6) and (7). The shots are taken at different times t " 0, 1, 5, 10. It consists in stopping (value 1) when the process reaches the plastic boundary and continuing otherwise (value 0). The terminal condition is f " 1 t|z|"1u , the right hand side is g " 0 and λ " 0.1.
accumulated plastic deformation (see Figure 8). Second, we could consider more general control problems (in this work, as a first step we have considered only a control acting linearly on the drift, see (8)). In particular, computing solutions for the problems (2) remains a challenging task.    (10): evolution of the value function in time vpz, y, tq for the accumulated plastic deformation. The shots are taken at different times t " 0, 1, 5, 10.
The terminal condition is f " 0, the right hand side is g " |y|1 t|z|"1u and λ " 0.1.  (10): evolution of the optimal control in timeâpz, y, tq for the accumulated plastic deformation. The shots are taken at different times t " 0, 1, 5, 10.
The terminal condition is f " 0, the right hand side is g " |y|1 t|z|"1u and λ " 0.1. . HJB problem as shown in (10): evolution of the optimal control in timeâpz, y, tq for the accumulated plastic deformation (view from above). The shots are taken at different times t " 0, 1, 5, 10. The terminal condition is f " 0, the right hand side is g " |y|1 t|z|"1u and λ " 0.1.  Figure 10. HJB problem as shown in (10): Left: one realization of a trajectory with the optimal control (full line) and without control (dashed line) for the same realization of the noise; right: corresponding evolution of the optimal control in timeα t "âpz t , y t , tq for the accumulated plastic deformation. The terminal condition is f " 0, the right hand side is g " |y|1 t|z|"1u and λ " 0.1.