Combined Plant and Control Co-design via Solutions of Hamilton-Jacobi-Bellman Equation Based on Physics-informed Learning
Kenjiro Nishimura Hikaru Hoshino Eiko Furutani






Combined Plant and Control Co-design via Solutions of Hamilton-Jacobi-Bellman Equation Based on Physics-informed Learning




    
Kenjiro Nishimura,   Hikaru Hoshino,   Eiko Furutani
Department of Electrical Materials and Engineering 

University of Hyogo

2167 Shosya, Himeji, Hyogo 671-2280, Japan 

er24h020@guh.u-hyogo.ac.jp, {hoshino, furutani}@eng.u-hyogo.ac.jp 













    
================================================================================================================================================================================================================================================================



§ ABSTRACT 
This paper addresses integrated design of engineering systems, where physical structure of the plant and controller design are optimized simultaneously. To cope with uncertainties due to noises acting on the dynamics and modeling errors, an Uncertain Control Co-design (UCCD) problem formulation is proposed. Existing UCCD methods usually rely on uncertainty propagation analyses using Monte Calro methods for open-loop solutions of optimal control, which suffer from stringent trade-offs among accuracy, time horizon, and computational time. The proposed method utilizes closed-loop solutions characterized by the Hamilton-Jacobi-Bellman equation, a Partial Differential Equation (PDE) defined on the state space. A solution algorithm for the proposed UCCD formulation is developed based on PDE solutions of Physics-informed Neural Networks (PINNs). Numerical examples of regulator design problems are provided, and it is shown that simultaneous update of PINN weights and the design parameters effectively works for solving UCCD problems. 



Control Co-design (CCD), optimal control, Physics-informed Neural Networks (PINNs), stochastic processes





§ INTRODUCTION
 


When designing an engineering system, one can adopt a sequential strategy where the physical system design is optimized first, followed by the controller design. 
For instance, consider the system-level design of a robotic manipulator.
The physical system design may optimize the geometric parameters of the links, and the control design may determine the joint torque time trajectories for specific tasks <cit.>.
Although the physical system design is often performed based only on static characteristics of the system, the resulting systems are suboptimal in most cases, and it is desirable to consider dynamic characteristics in the physical system design in cooperation with the controller design.  
This view is called as Control Co-design (CCD) approach <cit.>, and many authors have shown its benefit in various applications including robotic manipulators <cit.>, quadruped robots <cit.>, flexible space structures <cit.>, electric motors <cit.>, offshore wind farms <cit.>, Field Programmable Gate Array (FPGA) circuits <cit.>, and civil structures <cit.>. 

One of the main challenges in CCD that remains to be addressed is to consider the impact of uncertainties coming from the noise acting on the control channels, unmodeled or neglected dynamics of the system, estimation errors in model parameters, and so on. 
All of these uncertainties may propagate through the dynamical system and transform the states into uncertain trajectories <cit.>.
These problems are termed as Uncertain CCD (UCCD) <cit.> or Robust CCD <cit.>, and several formulations and solution methods have been proposed. 

In <cit.>, it is proposed to optimize a metric that represents the sensitivity of the trajectory to perturbation.
This can reduce the sensitivity of the trajectory for a specific uncertainty relying on custom-made cost formulations, but it increases the complexity of the problem and the cost function has to be carefully designed.  
Another approach is stochastic programming, in which the optimal trajectory is found for a set of perturbed scenarios <cit.>.
In <cit.>, a bi-level optimization scheme, where trajectories are optimized by the Differential Dynamic Programming (DDP) algorithm in an inner loop, and hardware is optimized in an outer loop with a genetic algorithm based on Covariance Matrix Adaptation Evolution Strategy (CMA-ES) <cit.>.
However, the above methods rely on Monte Calro simulation for analyzing uncertainty propagation through dynamical systems, and are subject to stringent trade-offs among accuracy, time horizon, and computational time <cit.>. 


In this paper, we propose a novel approach to UCCD problems, where optimality of the trajectories are characterized by the Hamilton-Jacobi-Bellman (HJB) equation in optimal control. 
In contrast to DDP based approach in <cit.>, where open-loop solutions of an optimal control problem are computed, the HJB equation characterizes closed-loop solutions in the form of a deterministic Partial Differential Equation (PDE), and uncertainty propagation can be naturally incorporated within the framework of stochastic control theory <cit.>. 
In our previous work <cit.>, we proposed a CCD solution method for deterministic dynamics with uncertainties in initial conditions based on Galerkin-approximations of the HJB equation <cit.>. 
This paper presents an extension of the CCD framework in <cit.> to deal with stochastic dynamics and uncertainties in model parameters. 
To the best of our knowledge, this is the first work to use the HJB equation for the purposes of solving UCCD problems. 
Furthermore, we present a UCCD solution method based on Physics-informed Neural Networks (PINNs) <cit.>. 
PINNs are widely used for solving various PDEs <cit.> and have a potential in dealing with a high-dimensional system <cit.>, whereas the HJB equation is notoriously difficult to solve by using standard numerical methods when the dimension becomes 5 or more <cit.>. 
This paper presents numerical examples of regulator design for a 2-dimensional nonlinear system and Linear Quadratic Regulator (LQR) problems of up to 10 dimensional systems to examine the effectiveness of the proposed method. 

The rest of this paper is organized as follows. 
In sec:formulation, the proposed UCCD problem formulation is presented. 
A solution method for the UCCD problem based PINNs is introduced in sec:method. 
Numerical results are provided in sec:simulation. 
Conclusions and future works are summarized in sec:conclusion. 




 §.§ Notation


Let ℝ be the set of real numbers, and ^n be the n-dimensional Euclidean space. 
For an open set A, A stands for the closure of A, and ∂ A for the boundary of A. 
For a scalar function ϕ, 
∂_x ϕ stands for the gradient of ϕ with respect to x, and ∂_x^2 ϕ for the Hessian matrix of ϕ. 
Let tr(M) be the trace of the matrix M. 



For random variables X and Y, let 𝔼[X] be the expectation of X, and 𝔼[X|Y=y] be the conditional expectation of X given Y=y. 
We use upper-case letters (e.g., Y) to denote random variables and lower-case letters (e.g., y) to denote their specific realizations. 




§ CO-DESIGN PROBLEM FORMULATION
 

In this section, we firstly introduce a basic co-design problem commonly studied in literature in sec:basic_formulation, and then present the proposed formulation in sec:proposed_formulation. 




 §.§ Deterministic Control Co-design Problem
 

There are two commonly studied strategies for CCD problems: simultaneous and nested <cit.>.
The simultaneous strategy optimizes both the plant and control variables in a same optimization formulation and is the most fundamental representation <cit.>. 
The nested strategy can be seen as a specific reorganization of the simultaneous formulation. An outer loop optimizes the design of the controlled plant, and an inner loop identifies the optimal control for each plant design tested by the outer loop. 
Here we briefly introduce a basic formulation in the simultaneous strategy. 
Consider the following dynamical system with an l-dimensional design parameter ρ=[ρ_1,  …,  ρ_l ]^⊤: 

    ẋ = f( x, u; ρ )

where x ∈ℝ^n stands for the n-dimensional state, ẋ for the time derivative of x, and u ∈ℝ^m for the m-dimensional control action.  
For this system, a standard optimal control problem can be considered with the following cost integral J_c: 

    J_c(ρ, u(·)) =  
     ∫_0^T  L ( x(t), u(t), ρ) dt 
     + M(x(T), ρ)

where T stands for the horizon length, and L and M for the Lagrange cost (running cost) and Mayer cost (terminal cost), respectively.  
The above notation shows that the cost integral J_c can be affected by the design parameter ρ through the dependence of L and M on ρ. 
Besides, J_c depends on ρ through the change in the dynamics (<ref>). 
By combining the cost J_c representing the control performance and an additional cost J_p about the choice of the parameter ρ of the plant, which represents, e.g., hardware materials costs or assembling costs, the co-design problem can be formulated as 

    min_ρ,  u(·)    w_p J_p(ρ) + w_c J_c(ρ, u(·))  
       s.t.  ẋ(t) = f(x(t), u(t); ρ),   ∀ t∈ [0, T],  
         x(0) = x_0

where w_p and w_c are weighting coefficients, and x_0 stands for the initial state of the optimal control problem. 
Although a more concise formulation can be obtained by including the term J_p in the Mayer term of the optimal control problem, the above formulation is commonly used to allow for more natural representations of CCD problems. 





 §.§ Proposed Formulation
 


Here we present the proposed formulation for UCCD problems. 
Let (Ω, ℱ, {ℱ_t}_t≥0, ) be a filtered probability space, and consider a control system with stochastic noise represented by ℱ_t-standard w-dimensional Brownian motion { W_t }_t ≥ 0 starting from W_0 = 0.
For an open set 𝕏⊂^n in the n-dimensional state space, the state X_t ∈𝕏 evolves according to the following Stochastic Differential Equation (SDE): 

    dX_t = { f(X_t;ρ)+g(X_t;ρ)U_t }dt + σ(X_t; ρ) dW_t,

where ρ stands for the design parameter of the hardware, and {U_t}_t≥0 with U_t ∈𝕌⊂ℝ^m is 
an {ℱ_t}_t≥0-adapted control process, which takes values in 𝕌. 
Throughout this paper, we assume sufficient regularity in the coefficients of the system (<ref>).
That is, the functions f, g and σ are chosen in a way such that the SDE (<ref>) admits a unique strong solution (see, e.g., Section IV.2 of <cit.>). 
The size of σ(X_t; ρ) is determined from the uncertainties in the disturbance, unmodeled dynamics, and prediction errors of the environmental variables. 
As we solve the equation in the domain 𝕏, define a stopping time τ as 

    τ = inf{ t  |   X_t ∉𝕏}.

Then, for each ρ, consider an optimal control problem to minimize the following cost functional: 

    J_c(ρ, x, U) =     [ ∫_0^τ L(X_s, U_s, ρ) e^-γ sds 
        + e^-γτ M(X_τ, ρ)  |  X_0 = x  ],

where L and M represent Lagrange and Meyer costs as mentioned in sec:basic_formulation, and γ stands for the discount factor. 
The control process U:={U_t}_t≥0 is chosen over a set 𝒰 of admissible control processes that have values in 𝕌 and are adapted to the filtration {ℱ_t}_t≥0. 
Then, by defining the value function V as 

    V(ρ, x) = inf_U∈𝒰 J_c (ρ, x, U),

the control performance is characterized as a function of the parameter ρ and the initial state x. 
From stochastic control theory <cit.>, it can be shown that the value function (<ref>) satisfies the following HJB equation:


    inf_u ∈𝕌{ℒ^u V(ρ, x) + L(x,u,ρ) -γ V(ρ, x) }=0,

in 𝕏, where ℒ^u is defined by 

    ℒ^u V(ρ,x) := 
           12Tr[ σ^⊤σ (∂_x^2 V) ] (ρ, x) 
        + { f(x;ρ)+g(x;ρ)u}^⊤∂_x V(ρ,x)

with the boundary condition 

    V(ρ, x) = M(ρ, x),

on ∂𝕏. 
By using the above, the UCCD problem proposed in this paper is formulated as

    min_ρ  J = w_p  J_p(ρ)  + w_c∫_𝕏ω(x) V(ρ,   x) dx  
       s.t.
       inf_u ∈𝕌{ℒ^u V(ρ, x) + L(x,u,ρ) -γ V(ρ, x) }=0,    x ∈𝕏,   
        V(ρ, x) = M(ρ, x),    x ∈∂𝕏

where ω: ℝ^n →ℝ is a weighting function to take an expectation of the value V(ρ, x) satisfying

    ∫_𝕏ω(x) dx = 1.

The problem (<ref>) is a partial differential equation constrained optimization problem, and a proper computational technique is needed to be solved. 



 In the above formulation, both the effects of the noise acting on the control channels and unmodeled/neglected dynamics are represented by the nonlinear function σ. 
 Also, uncertainties in initial conditions are represented by the distribution ω. Furthermore, if there are uncertain model parameters, denoted by ϕ, the system dynamics can be augmented as 
 
    [ dX_t;   dϕ ]= 
          [ f(X_t;ρ, ϕ)+g(X_t; ρ, ϕ) U_t;                            0 ]dt
       
        +[ σ(X_t; ρ);         0 ]dW_t.

 Thus, by applying the formulation (<ref>) adopted for the augmented system (<ref>), uncertainties in model parameters can also be captured by the distribution ω. 
 In conclusion, the proposed formulation can capture all the uncertainties due to the noise acting on the control channels, unmodeled/neglected dynamics, initial conditions, and model parameters. 






§ SOLUTION ALGORITHM
 







This section presents a solution method for the proposed formulation for UCCD problems. 
To deal with the HJB equation as a PDE constraint in the optimization problem (<ref>), we use the PINN framework, which is able to solve PDEs by exploiting machine learning techniques <cit.>.  
fig:method shows a schematic overview of the proposed UCCD method based on PINN. 
The PINN takes the pair (ρ, x) of the design parameter ρ and the state x (as well as the uncertain model parameter ϕ stated in Remark 1 if the augmented representation is used), and outputs the prediction V of the value V(ρ, x), as well as its derivatives ∂V/∂ρ and ∂V/∂ x computed by automatic differentiation <cit.>. 
By assuming that the PINN is parameterized by θ, the loss function L_PINN for the learning is defined as 
    L_PINN(θ, 𝒮_h, 𝒮_b) = μ_h L_HJB(θ, 𝒮_h) + μ_b L_bdry(θ,𝒮_b)

where L_HJB and L_bdry are loss terms for the HJB equation (<ref>) and the boundary condition (<ref>), respectively, and μ_h and μ_b are the weighting coefficients.
The HJB loss term L_HJB is a function of θ and a set of N_h random samples with 𝒮_h = { (ρ_i, x_i)  |  i∈{1,… N_h}, x_i ∈𝕏}, and  given by 

    L_HJB(θ, 𝒮_h) = 1N_p∑_i=1^N_pF(ρ_i, x_i, θ)   ^2

with 
 
    F(ρ, x, θ)    = inf_u ∈𝕌{ℒ^u V(ρ, x, θ) + L(x,u,ρ) -γV(ρ, x) }.

Here, the optimal control u needs to be addressed, and in this paper, we assume a specific form of the cost function. 

    The Lagrange cost term L in the cost functional J_c in (<ref>) takes the following form  
    
    L(x, u) = L̂(x) + u^⊤𝖱 u,

    where 𝖱 is a positive-definite matrix. 

With this assumption, we have an explicit optimal control as 

    u^∗(x) = -12𝖱^-1 g(x)^⊤∂_x V.


The boundary loss term L_bdry in (<ref>) is computed from a set of N_b random samples with 𝒮_b = { (ρ_j, x_j)  |  j∈{1,… N_b}, x_j ∈∂𝕏}, and given by 

    L_bdry(θ,𝒮_b) = 1N_b∑_j=1^N_bV(ρ_j, x_j, θ) - M(ρ_j, x_j) ^2.



To solve the UCCD problem (<ref>), we propose to simultaneously update the weights θ of the PINN and the design parameter ρ. 
The proposed solution method is presented in Algorithm <ref>. 
At each epoch, the sets 𝒮_h = { (ρ_i, x_i )}_i=1^N_h and 𝒮_b = { (ρ_j, x_j )}_j=1^N_b are randomly sampled with x_i ∈𝕏 and x_j ∈∂𝕏. 
The samples ρ_i and ρ_j are drawn by 

    ρ_i=ρ+ϵ_i,   ρ_j = ρ + ϵ_j,   ϵ_i,ϵ_j ∼𝒩(0, σ_n^2)

where 𝒩(0, σ_n^2) stands for the Gaussian distribution with the zero mean and the standard deviation of σ_n. 
The noise term ϵ is added for exploration of ρ, and its effect is discussed in sec:simulation. 
The minimization of the loss L_PINN imposes the constraints in the UCCD problem (<ref>). 
The design parameter ρ is then updated to minimize the objective function J.  
To this end, a set 𝒮_r of N_r random samples are drawn as 𝒮_r = { x_k  |,  k∈{ 1, …, N_r}, x_k ∈𝕏}, and the following loss L_r is considered:

    L_r(ρ, 𝒮_r) =  w_p  J_p(ρ)  + w_c∑_k=1^N_rV(ρ,   x_k)

Note that the samples x_k needs to be drawn from the distribution ω in (<ref>) to be unbiased. 
The parameter ρ is updated by performing a gradient step on the loss L_r once in N_up epochs. 










§ NUMERICAL EXAMPLES
 

This section provides numerical examples of CCD problems for regulator designs.
A 2-dimensional nonlinear deterministic dynamical system is treated in sec:planer, and LQR problems for up to 10-dimensional systems are treated in sec:stochasticLQR. 
The algorithm was implemented by a deep learning framework PyTorch <cit.>, and our implementation is available at <https://github.com/er24h020/SCIS_ISIS_2024>. 




 §.§ Deterministic Nonlinear Planar System
 

We consider a CCD problem for a planer nonlinear dynamical system with uncertainty in initial states, for which a numerical analysis has been performed in our previous work <cit.> based on a Galerkin approximation-based CCD method.  
Consider the following deterministic dynamics:

    [ ẋ_1; ẋ_2 ]_ẋ 
     = 
     [ -x_1^3 -x_2;    x_1 +x_2 ]_f(x; ρ)
     + [ 0; ρ ]_g(x; ρ)
     u

where x = [x_1, x_2]^⊤∈^2 stands for the state, u ∈ for the input, and ρ∈ℝ for the design parameter. 
The objective function for the optimal control is given by 

    J_c(ρ, x) = ∫_0^∞ p x(t)^⊤ x(t) + q u^2(t) dt,

with p=q=1.
With this cost function, we have the following HJB equation:

    ∂_x V^⊤ (f + g u)  +  p x^⊤ x+ q u^2 =0

with u= -ρ∂_x_2V/(2q). 
The boundary condition is given as V(ρ, 0) = 0. 
Then, consider the CCD problem with the following objective function:

    J = w_p |ρ| + w_c∫_𝕏ω(x) V(x) dx

where 𝕏 =  { (x_1, x_2) ∈ℝ^2  |  |x_1|≤ 1, |x_2|≤ 1  }, w_p=1,  w_c=4, and ω(x)≡ 1/4. 
With this problem setting, the control performance is improved as ρ increases, but a penalty is imposed on the increase in ρ by the first term of the objective function. 

The objective function J was minimized by Algorithm <ref>. 
For the function approximator V of the value function, we used a neural network with 3 hidden layers with 64 units per layer, and the hyperbolic tangent () as the activation function. 
The soft plus function was used at the output. 
For updating the neural network's weights,  optimizer was used with the learning rate of 3E-3. 
The weights in the loss (<ref>) was set as μ_h=μ_b=1, and the number of samples was N_p=1000 for the inside of 𝕏 and N_b=100 for the boundary. 
The design parameter ρ was initially set to ρ = 1, and updated by using  optimizer with the learning rate of 2e-2 at every 1 or 500 epochs (N_up=1 or 500) with N_r=1000. 







fig:2d_results shows results of the proposed CCD algorithm applied to the planer example.
We compare the cases where ρ is sampled with and without adding the noises ϵ_i and ϵ_j at the lines 4 and 6 in Algorithm <ref>. 
For the additive noise, we used a uniform distribution with the amplitude of 0.1. 
In fig:rho_nolta, cyan and blue lines show the changes in the design parameter ρ when it is updated at each epoch (N_up=1). The purple line show the results with N_up=500 and adding noise. 
The solid curves correspond to the mean of ten repeated experiments, and the shaded region shows their standard deviations.
It can be seen that the standard deviation of ρ is significantly reduced by adding the noise, and the parameter ρ converges to the optimal solution ρ = 2.1, which is calculated in <cit.>, only when ρ is explored by adding the noise. 
As shown in fig:loss_nolta,  the PINN loss L_PINN is kept small, and it can be confirmed that the simultaneous update of ρ does not impede the learning process of the PINN. 


    



 §.§ Stochastic LQR problem
 

Next we present a stochastic LQR example based on <cit.>. 
Consider the following d-dimensional controlled stochastic process given by 

    dX_t = ρ U_t dt + √(2)dW_t

where X_t ∈𝕏⊂^d, U_t ∈^d, W_t ∈^d, and ρ∈. 
For the domain 𝕏, consider d-dimensional sphere of radius R: 

    𝕏 = { x ∈^d  |   x  < R }.

The cost functional for the optimal control is given by 

    J_c(ρ, x) =  [    ∫_0^τ ( p  X_t ^2 + q   U_t ^2 -2k d) e^-γ sds  
        +  e^-γτ k R^2 ]

where τ stands for the exit time of the domain 𝕏, and the constant k is given by 
  k = (√(q^2 γ^2 + 4pqρ^2) -γ q )/( 2 ρ^2 ). 
With this setting, the value function V satisfies the following HJB equation:

    ∂_x^2 V(ρ, x) + inf_u ∈^d  (ρ u^⊤∂_x V(ρ,x) + qu^2 ) 
        + px^2 -2kd -γ V(ρ, x) = 0.

with the boundary condition V(ρ, x) = k R^2. 
It is known that the PDE has the exact solution as a quadratic function V(x) =k  x^2, and the optimal control is given as u^∗(x)= -ρ∂_x V(x)/(2q) = -kρ x /q. 
Then, the objective function for the UCCD problem is given by 

    J = w_p |ρ| + w_c∫_𝕏ω(x) V(x) dx

where w_p= w_c=1, and ω(x)≡ 1/|𝕏|. 







fig:lqr_results shows results for d=5 (shown by blue) and 10 (shown by purple). The other parameters are given as p=q=γ=1, and R=2 and 4 for d=5 and 10, respectively. 
The shape of the neural network is the same with the above example, and the learning rates are 1e-4 and 1e-2 for θ  and ρ, respectively, for both of d=5 and 10.  
The weights in the loss (<ref>) are μ_p=μ_b=1. 
The numbers of samples are N_p = N_b = N_r = 1000 for d=5, and N_p = N_b = N_r = 100,000 for d=10. 
fig:lqr_rho shows the mean (solid lines) and standard deviation (shaded area) of the design parameter ρ in 10 repeated experiments, and it can be seen that ρ converges at around 120,000 epochs.
fig:lqr_loss shows the changes in the loss function L_PINN, and it is kept low for d=5 as shown by the blue line. 
For d=10, purple shows the mean value of the ten experiments, and cyan shows the results when the loss function is the smallest among these experiments. 
Although a relatively large variation was observed in the loss, the design parameter ρ converged to similar values in these 10 experiments. 









§ CONCLUSIONS
 

This paper proposed a novel UCCD problem formulation to cope with  uncertainties coming from noises acting on the dynamics and modeling errors.
The proposed method utilizes closed-loop solutions of an optimal control problem characterized by the Hamilton-Jacobi-Bellman equation as a PDE constraint in the UCCD problem.  
 A solution algorithm is developed based on Physics-informed Neural Networks (PINNs), and numerical examples show that simultaneous update of PINN weights and the design parameters effectively works for solving UCCD problems.

Future directions of this work include extension of the proposed algorithm to address problems where the optimal control is obtained only in an implicit form. 
A possible approach for this task is to use a reinforcement learning framework to obtain an optimal policy. 
Another future direction is to used the proposed method in practical applications in e.g., robotics and energy systems as mentioned in sec:introduction. 





§ ACKNOWLEDGMENT


This work was supported in part by JST, ACT-X Grant Number JPMJAX210M, Japan, and the Kansai Research Foundation for Technology Promotion, Japan. 




IEEEtran












