Sub-optimal Tracking in Switched Systems with Controlled Subsystems and Fixed-mode Sequence using Approximate Dynamic Programming
Tohid Sardarmehni Xingyong Song
eess.SY cs.SY

[NO \title GIVEN]
    [NO \author GIVEN]
    December 30, 2023
======================





 
Optimal tracking in switched systems with controlled subsystem and Discrete-time (DT) dynamics is investigated. A feedback control policy is generated such that a) the system tracks the desired reference signal, and b) the optimal switching time instants are sought. For finding the optimal solution, approximate dynamic programming is used. Simulation results are provided to illustrate the effectiveness of the solution.









§ INTRODUCTION
 

In this study, switched systems are referred to as systems comprised of several subsystems/modes such that at each time instant only one subsystem is active. If the subsystems include continuous control signals, the subsystems are called controlled subsystems. In case no control input exists in the subsystems, the subsystems are called autonomous subsystems. Control of switched systems with controlled subsystem is a challenging problem as the controller needs to define the switching time and the continuous control in the subsystems. Also, if the sequence of active subsystems. i.e., mode sequence, is known a priori, the mode sequence is called a fixed mode sequence. In a fixed mode sequence, the controller assigns the switching times and does not assign the active modes. 

On the other hand, if the mode sequence is free, then the controller needs to decide which mode should be activated at each time instant. Due to the discontinuous nature of the switched systems, control of these systems are challenging. Control of switched systems addresses many critical problems in automotive engineering <cit.>, power engineering <cit.>, and many other engineering fields.  


From the mathematical point of view, solutions of the underlying Hamilton-Jacobi-Bellman (HJB) equation provide the necessary and sufficient condition for optimality in optimal control problems <cit.>. However,  solving the HJB equation explicitly is difficult and most cases impossible. 

Dynamic Programming (DP) can solve optimal control problems through discretizing the state and control domain and finding the optimal value functions backward in time. However, as the order of the system increases, rapid access to memory becomes prohibitive in DP which is known as curse of dimensionality. To remedy the curse of dimensionality, one solution is to seek a near-optimal control solution instead of the exact optimal control which is provided by Approximate Dynamic Programming (ADP). In general, ADP methods use function approximators to approximate the value functions and they use iterative methods to tune the parameters of these function approximators.

ADP solutions for optimal control of switched systems with free mode sequence were investigated in <cit.>. As for the fixed mode sequence, a transformation was introduced in <cit.> to treat the switching times as parameters. This transformation was adapted in <cit.> and the mode sequence was incorporated in the value functions to solve a regularization problem. 



To solve the optimal tracking problem in systems with conventional dynamics, some ADP solutions were introduced in <cit.>. In <cit.>, an ADP solution based on Single Network Adaptive Critic (SNAC) structure was introduced which approximated the optimal costates and used them to provide the optimal policy. Also, another ADP solution was developed in <cit.> which approximated the value functions. For this method, a change of variables was performed in the system and a new state vector was introduced which included the tracking error dynamics and the dynamics of the reference signal. Furthermore, in <cit.> model free tracking methods were investigated. As for the switching dynamics, optimal tracking was studied in <cit.> in switched systems with autonomous subsystems and free mode sequence, and in <cit.> in switched systems with controlled subsystems and free mode sequence.    

In this paper, optimal tracking in switched systems with controlled subsystems and fixed mode sequence is studied. The basic idea in this study is using the transformation introduced in <cit.> to parametrize the switching instants and then use the idea introduced in <cit.> to include the mode sequence in the value functions and costates. Hence, the present study combines the results of <cit.> to develop a solution for optimal tracking in switched systems with fixed mode sequence. To solve the optimal tracking problem, a single network is used to capture the optimal costates which are parametrized with respect to the switching instants. Once the optimal costates are known, some recommendations are given to find the optimal switching times from the costates at an initial condition.
  

The rest of the manuscript is organized as follows. In section <ref>, optimal control problem formulation and some assumptions are presented. In section <ref>, the proposed solution based on the SNAC is introduced. 

Simulation results are discussed in section <ref>, and section <ref> concludes the paper.  





§ PROBLEM FORMULATION
 

The dynamics of the switched systems can be shown as  


    ẋ(t)= f̅_v(x(t)) + g̅_v( x(t) ) u(t),   v∈𝒱 = {1,2,…, M}, x(0) = x_0


where x∈ℝ^n is the state vector, u∈ℝ^m is the input, and t denotes the time. The Lipschitz continuous functions f̅_v: ℝ^n →ℝ^n and g̅_v: ℝ^n →ℝ^n× m denote the dynamics of the subsystems. The sub-index v portrays the active mode which can be selected from the set of all available modes, 𝒱, in the system. It is further assumed that f̅_v(0) = 0, for all modes v∈𝒱. The inclusion of continuous control, i.e., u(.), in (<ref>) shows that the subsystems are controlled subsystems.  

Assuming the sequence of active modes is known, it is desired to find the continuous control policy u(.), and the switching times, such that a performance index presented as 


    J   (x_0 )= (x(t_f)- r(t_f))^T S ( x(t_f) - r(t_f) ) +  ∫_t_0^t_f1/2((x(t)-r(t))^T Q̅(x(t)-r(t))+ u(t)^T R̅ u(t) ) dt


is minimized. In (<ref>), t_0 is the initial time, and t_f is the final time. r ∈ℝ^n is the reference signal where ṙ(t) = f̅_r_v(t) is the dynamics of the reference signal and the sub-index r_v shows the active reference dynamics. S ∈ℝ^n× n is a positive semi-definite matrix for penalizing the terminal cost, Q̅∈ℝ^n× n is the state penalizing matrix which is assumed to be positive semi-definite, and R̅∈ℝ^m × m is a positive definite control penalizing matrix. 



Using Euler integration method, by choosing a small sample time δ t>0, one can discretize the dynamics (<ref>) as


    x_k+1 = f_v(x_k) + g_v(x_k) u_k


where the non-negative integer k is the discrete time index. For notational simplicity, the discrete time index is shown as a sub-index, i.e., x_k ≡ x(k). Also, f_v(x_k) = x_k + f̅_v(x_k) δ t, and g_v(x_k) = g̅_v(x_k) δ t. With a similar procedure, one can discretize the cost function (<ref>) as 


    J   (x_0)= (x_N - r_N)^T S ( x_N - r_N )  +  ∑_k = 0 ^N1/2((x_k-r_k)^T Q (x_k-r_k)+ u_k^T R u_k)


In (<ref>), N = t_f - t_0/δ t, Q = Q̅δ t, and R = R̅δ t. Based on (<ref>), one can define the cost-to-go as the cost of going from discrete time index k to N as 


    J   (x_k)= (x_N- r_N)^T S ( x_N - r_N ) +  ∑_k̅ = k ^N1/2( (x_k̅-r_k̅)^T Q (x_k̅-r_k̅)+ u_k̅^T R u_k̅)


Before going forward, the following definition and assumption are needed. 

 
A control policy is called admissible if it stabilizes the system presented in (<ref>) in a selected compact region of interest Ω⊂ℝ^n, which includes the origin. Also, ∀ x_0 ∈Ω, J(x_0) should be finite if the state is propagated along that policy.


 
Given the mode sequence, there is at least one admissible policy for the system. 



Considering Assumption <ref>, one can define the value function as 

    V   (x_k, k)≡ V_k(x_k)= min_u(.)((x_N- r_N)^T S (x_N - r_N) + 1/2∑_k̅ = k ^N( (x_k̅-r_k̅)^T Q (x_k̅-r_k̅)+ u_k̅^T R u_k̅) )


Considering time step k to k+1, one has 


    V_k(x_k)= min_u(.)((x_N- r_N)^T S ( x_N - r_N ) + 1/2(x_k-r_k)^T Q (x_k-r_k)+ 1/2u_k^T R u_k + 1/2∑_k̅ = k+1 ^N (x_k̅-r_k̅)^T Q (x_k̅-r_k̅)+ u_k̅^T R u_k̅)


After some algebraic manipulations, one can rewrite (<ref>) as


    V_k(x_k)= min_u(.)(   1/2(x_k-r_k)^T Q (x_k-r_k)+ 1/2u_k^T R u_k + V_k+1(x_k+1))
 

Equation (<ref>) simply means minimum cost of going from time k to N equals to cost of going from time k to k+1 plus minimum cost of going from time k+1 to N. This is in fact the Bellman equation of optimality <cit.>. Based on (<ref>), one can define the optimal policy as 


    u_k(x_k)= _u(.)(   1/2(x_k-r_k)^T Q (x_k-r_k)+ 1/2u_k^T R u_k + V_k+1(x_k+1))
 









 §.§ Including the Mode Sequence


In <cit.>, a transformation is introduced to include the switching times as parameters in the optimal control formulation. For ease of presentation, a case with two subsystems and only one switching is considered. Let the switching happen at t = t_1∈ (t_0,t_f). Also, let the mode sequence to be {mode 1, mode 2}[When t< t_1 mode 1 is active and when t≥ t_1 mode 2 is active. ]. To make the switching time instant an independent parameter, let <cit.> 


    t =
    
    t_0 + (t_1 - t_0) t̂    if 0 ≤t̂< 1
    
    t_1 + (t_f - t_1) (t̂ - 1)     if 1≤t̂≤ 2


From the transformation introduced in (<ref>), one notices that t ∈ [t_0, t_f] and t̂∈ [0,2]. The merit of the transformation is that the switching time t_1 can be any point in t∈ [t_0, t_f]. However, in the transformed time, i.e., t̂∈ [0,2], switching only happens at t̂ = 1. For t̂≥ 1, mode 2 is active and for t̂ < 1, mode 1 is active. Based on the introduced transformation in (<ref>), one has



    x'(t̂) = dx/d t̂ =  dx/dtdt/dt̂


Since the mode sequence is known, (<ref>) becomes


    x'(t̂)  =
    (f̅_1(x(t̂)) + g̅_1(x(t̂))u(t̂) ) (t_1 - t_0)     if 0 ≤t̂< 1
    (f̅_2(x(t̂)) + g̅_2(x(t̂))u(t̂) ) (t_f - t_1)     if 1≤t̂≤ 2

Also, the cost function in (<ref>) can be written as[Since the mode sequence is known, one can consider the integral from t_0 to t_1 with the first mode, and from t_1 to t_f with the second mode. ]


    J(x_0)  =   (x(2) - r(2))^T S (x(2) - r(2)) + ∫_0^11/2( (x(t̂) - r(t̂))^T Q̅ (t_1 - t_0) (x(t̂)- r(t̂)) + u(t̂)^T R̅(t_1 - t_0) u(t̂) )dt̂
        + ∫_1^21/2( ( x(t̂)- r(t̂))^T Q̅ (t_f - t_1) (x(t̂)- r(t̂)) + u(t̂)^T R̅(t_f - t_1) u(t̂) )dt̂


An important observation in (<ref>) is that the transformed cost function is not only a function of x_0, but also it is a function of the switching time, i.e., t_1. Hence, J(x_0) = J(t_1, x_0). With a similar procedure used before, by choosing a small sampling time δt̂ one can discretize (<ref>) and (<ref>) as


    x_k̂+1 =
    
    f_1(x_k̂)  + g_1(x_k̂) u_k̂    if 0 ≤k̂< N'/2
    
    f_2(x_k̂)  + g_2(x_k̂) u_k̂    if N'/2≤k̂≤ N'


where f_1(x_k̂ ) =  x_k̂ +  f̅_1(x_k̂) (t_1 - t_0) δt̂, g_1(x_k̂ ) =  g̅_1(x_k̂ ) (t_1 - t_0) δt̂, f_2(x_k̂ ) =  x_k̂ +   f̅_2(x_k̂) (t_f - t_1) δt̂, and g_2(x_k̂ ) = g̅_2(x_k̂) (t_f - t_1) δt̂. In (<ref>), k̂∈ [0 ,N'] is the discrete time index where N' = number of switching + 1/δt̂ <cit.>. 

With the transformed dynamics and the cost function, with similar procedure used in the previous section one can define cost-to-go as

 

    V_k̂(t_1, x_k̂) =  
    𝐐_1 + 𝐑_1 + V_k̂+1(t_1, x_k̂+1 )    if 0 ≤k̂< N'/2
    𝐐_2 + 𝐑_2 + V_k̂+1(t_1, x_k̂+1 )    if N'/2≤k̂≤ N'

where

    𝐐_1    = 1/2( x_k̂- r_k̂ )^T Q̅ (t_1 - t_0) δt̂ ( x_k̂- r_k̂)
    𝐑_1    = 1/2u_k̂^T R̅(t_1 - t_0)δt̂ u_k̂
    𝐐_2    = 1/2( x_k̂- r_k̂ )^T Q̅ (t_f - t_1) δt̂ ( x_k̂- r_k̂)
    𝐑_2    = 1/2u_k̂^T R̅(t_f - t_1)δt̂ u_k̂


As one can see, the cost-to-go in (<ref>) is a function of current time k̂, current state x_k̂, and the switching time t_1. Similarly, one can define the costate as[By definition, costate is the gradient of value function, i.e., λ(x) = ∂ V(x)/∂ x. Taking the gradient of (<ref>) leads (<ref>).]


    λ_k̂(t_1, x_k̂) =  
    𝒬̅_1 + ∂ x_k̂+1/∂ x_k̂λ_k̂ + 1(t_1, x_k̂+1)    if 0 ≤k̂< N'/2
    𝒬̅_2 + ∂ x_k̂+1/∂ x_k̂λ_k̂+1(t_1, x_k̂+1)    if N'/2≤k̂≤ N'


where 


    𝒬̅_1    = Q̅ (t_1 - t_0) δt̂( x_k̂- r_k̂)
    𝒬̅_2    = Q̅ (t_f - t_1) δt̂( x_k̂- r_k̂)







§ SINGLE NETWORK ADAPTIVE CRITIC (SNAC)
 

The application of SNAC was introduced for tracking in systems with conventional dynamics <cit.>. This idea is adapted in this section to perform tracking in a switched system. To introduce the concept, consider the costate as in (<ref>). 

The idea here is training neural networks to approximate λ_k̂+1(t_1, x_k̂+1) from (t_1,x_k̂). Based on Weierstrass Approximation Theorem <cit.>, linear-in-parameter neural networks with polynomial basis functions can uniformly <cit.> approximate continuous functions to a desired degree of precision in a compact set. In order to use Weierstrass Approximation Theorem, the following assumption is required. 

 
The value functions are continuously differentiable.


Through Assumption <ref>, one can use linear-in-parameter neural networks to approximate the value functions, and the costates. Consider the exact costate at discrete time index k̂ as 


    λ_k̂+1(t_1, x_k̂+1) = W_k̂^*^T ϕ(t_1, x_k̂) + ε_k̂^*(t_1, x_k̂)
 

where W_k^* ∈ℝ^m_λ× n is a weight vector and ϕ: ℝ×ℝ^n →ℝ^m_λ is a vector of linearly independent polynomial basis functions (neurons). The number of neurons is denoted by positive integer m_λ. In (<ref>), the dependence of the parameters/functions to discrete time index is shown with a sub-index k̂. Hence the approximate costates can be calculated as 


    λ_k̂+1(t_1, x_k̂+1) = W_k̂^T ϕ(t_1, x_k̂)


where W_k̂∈ℝ^m_λ× n is a tunable weight vector. The W_k̂ is tuned through the training process. Once the costates are known, one finds the optimal policy as 


    u_k̂(t_1, x_k̂) =  
    
    -R_1^-1g_1^T(x_k̂) λ_k̂+1(t_1, x_k̂+1)    if 0 ≤k̂< N'/2
    
    -R_2^-1g_2^T(x_k̂) λ_k̂+1(t_1, x_k̂+1)    if N'/2≤k̂≤ N'


In (<ref>), R_1= R δt̂(t_1 - t_0) and R_2 = R δt̂(t_f - t_1). For training, one can go backward in time and find the costates and save them for online control. This process is summarized in Algorithm <ref>. 






The convergence of the inner loop in step 3 of Algorithm <ref> was studied in theorem 1 of <cit.> for systems with conventional dynamics. The proof can be modified to include tracking in switched systems. 





Once the training is concluded, one needs to find the optimal switching times from the costates for a selected initial condition x_0 ∈Ω. Three methods are suggested below to find the optimal switching times from the optimal costates.  



  * – Method 1: propagating the states analytically along the optimal policy by treating switching time as a parameter and finding the optimal cost-to-go from the cost function. Once done, one can use constrained minimization methods to find switching times. 


  * – Method 2: integrating the costate analytically to find the value function. 
Similar to finding the velocity field from potential flow in fluid mechanics, one can integrate the costates analytically to find the value functions. The convenient feature of this method is that the analytical solutions provide the optimal value function ∀ x_0 ∈Ω. In order word, one does not need to integrate again when the initial condition is changed. This is unlike Method 1 that propagation should be done for each initial condition separately.  However, as the order of system increases, this method becomes very complicated. In other words, this method is only suitable for systems with low order dynamics.  


  * – Method 3: propagating the states along all possible switching times and find the optimal cost to go for all possible switching time. Once done, choose the switching times which lead to the minimum value function. Method 3 is similar to the forward dynamic programming method and when the number of switching increases, performing this method might become very time-consuming. 










§ NUMERICAL SIMULATION
 

In this section some simulation results are provided to evaluate the effectiveness of the solutions discussed in this paper. Consider a system with two modes. For the first mode, a benchmark system, Van der Pol oscillator was selected. The dynamics of this mode can be shown as 


    ẋ_1(t)    = x_2(t)
    ẋ_2(t)    = (1 - x^2_1(t) )  x_2(t) - x_1(t) + u(t)
 

For the second mode, a linear subsystem was selected as


    ẋ_1(t)    = x_2(t)
    ẋ_2(t)    = 2x_1 - x_2 + u(t)
 

The reference signal is chosen as 


    ṙ_1(t)    = sin(π t)
    ṙ_2(t)    = πcos(π t)
  




For the cost function, S = (10^5 , 10^5), Q̅ = (10^5, 10^7) were selected, where (a,b) denotes a diagonal matrix with a and b on the main diagonal and 0 elsewhere. Also, R̅ = 1/δt̂ was selected as the control penalizing term. The basis functions were selected as all possible combination of t_1, x_1, x_2 up to the power 3 without repetition. 


For training, 1000 random samples were generated in Ω = { (t_1, x_1, x_2) | t_1∈ (0,3), x_1∈ [-4,4],  x_2 ∈ [-4,4] }. For descritization, sampling time was selected as δt̂ = 0.001 (sec). Using approach 1, The training process concluded in 22.6177 seconds using an office desktop with 16 GB of RAM and Intel(R) Core(TM) i7-3770 CPU @ 3.4 GHz. The history of the weights of the neural networks to approximate the costates is shown in Fig. <ref>. As one can see from Fig <ref>, the history of the weights shows a jump at t̂ =1 which is the switching time.  

 


Once the training concluded, the trained costates were used to find the optimal switching time. For this example, Method 2 in Remark <ref>, i.e., integration of the costate analytically, was used to find the value function and then minimized the value function with nonlinear programming methods. Through integration, one can find the analytical value function as V(t_1,x_1, x_2) =67377 t_1^3 x_1 + 568.7 t_1^3 x_2 + 49288 t_1^2 x_1^2 + 849 t_1^2x_1x_2 - 2.322× 10^5t_1^2x_1 + 57.67t_1^2x_2^2 - 2077t_1^2x_2 + 4210t_1x_1^3 + 162.2t_1x_1^2x_2 - 92177t_1x_1^2 + 1.432t_1x_1x_2^2 - 1691t_1x_1x_2 + 7.968× 10^5t_1x_1 + 0.3332t_1x_2^3 - 2411t_1x_2^2 + 3337t_1x_2 + 2413x_1^4 + 97.95x_1^3x_2 - 4013x_1^3 - 434.8x_1^2x_2^2 - 182.2x_1^2x_2 + 1.44× 10^5x_1^2 - 2.424x_1x_2^3 + 2.186x_1x_2^2 + 1789x_1x_2 - 1.961× 10^6x_1 - 0.1953x_2^4 - 0.08602x_2^3 + 50288x_2^2 + 77666x_2. Substituting the initial conditions for x_1 = 1 and x_2 = -0.5, one can find the value function as a function of the switching time as V(t_1) = 67088t_1^3 - 1.823 × 10^5 t_1^2 + 7.073 × 10^5 t_1 - 1.846 × 10^6. Using nonlinear programming, the best switching time is sought as t_1 = 2.654 (sec). 

The history of the states in the transformed dynamics is shown in Fig. <ref>. As one can see, the controller had a very good performance in forcing x_2(.) to track r_2(.) as in the cost function it was emphasized.  

 






§ CONCLUSION
 
An approximate solution for optimal control of switched systems with fixed mode sequence and controlled subsystems was presented. The method includes two levels of control. In the upper level, optimal switching times were sought. In the lower level, continuous control for each mode was generated in a feedback form. To find the continuous control, a single network adaptive critic was used to find the optimal costates while treating the switching times as parameters. Simulation results confirmed the effectiveness of the solution.  




IEEEtran
