Online Optimal Control with Linear Dynamics and Predictions: Algorithms and Regret Analysis
Yingying Li Xin Chen Na Li


Fast uncertainty quantification of reservoir simulation with variational U-Net
    Larry Jin,   Hannah Lu,  Gege Wen

Energy Resources Engineering

Stanford University

{zjin, hannahlu, gegewen}@stanford.edu
    December 30, 2023
================================================================================================================================




	This paper studies the online optimal control problem with time-varying convex stage costs for a time-invariant linear dynamical system, where a finite lookahead window of accurate predictions of the stage costs are available at each time. We design online algorithms, Receding Horizon Gradient-based Control (RHGC), that utilize the predictions through finite steps of gradient computations. We study the algorithm performance measured by dynamic regret: the online performance minus the optimal performance in hindsight. It is shown that the dynamic regret of RHGC decays exponentially with the size of the lookahead window. In addition, we provide a fundamental limit of the dynamic regret for any online algorithms by considering linear quadratic tracking problems. The regret upper bound of one RHGC method almost reaches the fundamental limit, demonstrating the effectiveness of the algorithm. Finally, we numerically test our algorithms for both linear and nonlinear systems to show the effectiveness and generality of our RHGC.
	 




§ INTRODUCTION


In this paper, we consider a N-horizon discrete-time sequential decision-making problem. At each time t=0, …,N-1, the decision maker observes a state x_t of a dynamical system, receives a W-step lookahead window of future cost functions of states and control actions, i.e. f_t(x)+g_t(u), …, f_t+W-1(x)+g_t+W-1(u), then decides the control input u_t which drives the system to a new state x_t+1 following some known dynamics. For simplicity, we consider a linear time-invariant (LTI) system x_t+1=Ax_t+Bu_t with (A,B) known in advance. The goal is to minimize the overall cost over the N time steps. This problem enjoys many applications in, e.g. data center management <cit.>,  robotics <cit.>, autonomous driving <cit.>, energy systems <cit.>,  manufacturing <cit.>.

Hence, there has been a growing interest on the problem, from both control and online optimization communities. 






















In the control community, studies on the above problem focus on economic model predictive control (EMPC), which is a variant of model predictive control (MPC) with a primary goal on optimizing economic costs <cit.>. 
Recent years have seen a lot of attention on the optimality performance analysis of EMPC, under both time-invariant costs <cit.> and time-varying costs <cit.>. However, most studies focus on asymptotic performance and there is still limited understanding on the non-asymptotic performance, especially under time-varying costs. Moreover, for computationally efficient algorithms, e.g. suboptimal MPC and inexact MPC  <cit.>, there is limited work on the optimality performance guarantee.



In online optimization, on the contrary, there are many papers on the non-asymptotic performance analysis, where the performance is usually measured by regret, e.g., static regrets<cit.>, dynamic regrets<cit.>, etc., but most work does not consider predictions and/or dynamical systems. Further, motivated by the applications with predictions, e.g. predictions of electricity prices in data center management problems <cit.>, there is a growing interest on  the effect of predictions on the online problems <cit.>.
However, though some papers consider switching costs which can be viewed as a simple and special dynamical model <cit.>, there is a lack of study on the general dynamical systems and on how predictions affect the online problem with dynamical systems. 





























In this paper, we propose novel gradient-based online control algorithms, receding horizon gradient-based control (RHGC), and provide nonasymptotic optimality guarantees by dynamic regrets. RHGC can be based on many gradient methods, e.g.  vanilla gradient descent, Nesterov's accelerated gradient, triple momentum, etc., <cit.>. Due to the space limit, this paper only presents  receding horizon gradient descent (RHGD) and  receding horizon triple momentum (RHTM). For the theoretical  analysis, we assume strongly convex and smooth cost functions, whereas applying  RHGC does not require these conditions. Specifically, we show that  the regret bounds of RHGD and RHTM decay exponentially with the prediction window's size W, demonstrating that our algorithms efficiently utilize the prediction. Besides, our regret bounds  decrease when the system is more “agile” in the sense of a controllability index <cit.>. Further, we provide a fundamental limit for any online control algorithms and show that the fundamental lower bound almost matches the regret upper bound of RHTM. This indicates that  RHTM achieves near-optimal performance at least in the worst case. We also provide some discussion on the classic linear quadratic tracking problems, a widely studied control problem in literature, to provide more insightful interpretations of our results. Finally, we numerically test our algorithms. In addition to linear systems, we also apply  RHGC to a nonlinear dynamical system: path tracking by a two-wheeled robot. Results show that  RHGC works effectively for nonlinear systems though RHGC is only presented and theoretical analyzed on LTI systems. 

Results in this paper are  built on a paper on online optimization with switching costs <cit.>. Compared with <cit.>, this paper studies online optimal control with general linear dynamics, which includes  <cit.> as a special case; and studies how the system controllability index affects the regrets.






   There has been some recent work  on online  optimal  control problems with time-varying costs <cit.> and/or time-varying disturbances <cit.>,  but most papers  focus on the no-prediction cases. As we show later in this paper, these  algorithms can be used in our RHGC methods as initialization oracles. Moreover, our regret analysis shows that RHGC can reduce the regret of these no-prediction online algorithms by a factor exponentially decaying with the prediction window's size.


Finally, we would like to mention another related line of work:  learning-based control <cit.>. In some sense, the results in this paper are orthogonal to that of the learning-based control, because the learning-based control usually considers a time-invariant environment but unknown dynamics,  and aims to learn system dynamics or optimal controllers by data; while this paper considers a time-varying scenario with known dynamics but changing objectives and studies decision making with limited predictions. It is an interesting future direction to combine the two lines of work for designing more applicable algorithms.























Notations.
Consider matrices A and B, A≥ B means A-B is positive semidefinite and [A,B] denotes a block matrix. The norm · refers to the L_2 norm for both vectors and matrices. Let x^i denote the ith entry of the vector. Consider a set ={k_1, …, k_m }, then x^=(x^k_1,…, x^k_m)^⊤, and A(, :) denotes the  rows of matrix A stacked together. Let I_m be an identity matrix in ^m × m. 






































§ PROBLEM FORMULATION AND PRELIMINARIES


















Consider a finite-horizon discrete-time optimal control problem with time-varying cost functions f_t(x_t)+g_t(u_t) and a linear time-invariant (LTI) dynamical system:

    min_𝐱,𝐮      J(  𝐱,   𝐮)= ∑_t=0^N-1[ f_t(x_t)+   g_t(u_t) ] + f_N(x_N)
    s.t.      x_t+1=  A   x_t +   B   u_t,    t≥ 0

where  x_t ∈^n, u_t ∈^m, 𝐱=(x^⊤_1, …, x^⊤_N)^⊤, 𝐮=(u^⊤_0,…, u^⊤_N-1)^⊤, x_0 is given,   f_N(x_N) is the terminal cost.[The results in this paper can be extended to cost c_t(x_t,u_t) with proper assumptions.]  To solve the optimal control problem  (<ref>),  all  cost functions from t=0 to t=N are needed. However, at each time t, usually only a finite lookahead window of cost functions are available and the decision maker needs to make an online decision u_t using the available information. 


In particular, we consider a simplified prediction model: at each time t, the decision maker obtains accurate predictions for the next W time steps, f_t,g_t, …, f_t+W-1, g_t+W-1, but no further prediction beyond these W steps, meaning that  f_t+W,g_t+W,… can even be adversarially generated. Though this prediction model may be too optimistic in the short term and over pessimistic in the long term, this model i) captures a commonly observed phenomenon in predictions that short-term predictions are usually much more accurate than the long-term predictions; ii) allows researchers to derive insights for the role of predictions and possibly to extend to more complicated cases <cit.>. 



The online optimal control problem is described as  follows: at each time step t=0,1,…, 

		
  * the agent observes state x_t and receives prediction f_t, g_t,  …,  f_t+W-1, g_t+W-1;
		
  * the agent decides and implements a control u_t
		and  suffers the  cost f_t(x_t)+   g_t(u_t);
		
  * the system evolves to the next state x_t+1=Ax_t + Bu_t.[We assume known A, B, no process noises, state feedback, and leave relaxing assumptions as future work.]



	An online control algorithm, denoted as , can be defined as a mapping from the prediction information and the history information to the control action, denoted by u_t(𝒜):
	
    u_t()=(   x_t(), …,   x_0(), {  f_s,   g_s}_s=0^t+W-1),   t≥ 0,

	where x_t() is the state generated by implementing  and x_0()=x_0 is given.
	
This paper evaluates the performance of online control algorithms by comparing against the optimal control cost J^* in hindsight, that is, J^* min{ J(  𝐱,  𝐮)|   x_t+1=Ax_t+Bu_t,  ∀  t ≥ 0 }.

In this paper, the performance of an online algorithm     is measured by [The optimality gap depends on initial state x_0 and {f_t,g_t}_t=0^N, but we omit them for  simplicity of notation.]

    Regret(𝒜):= J(𝒜)-J^* = J(𝐱(), 𝐮())-J^*,

which is sometimes called as dynamic regret <cit.> or competitive difference <cit.>. 



Another popular regret notion is the static regret, which compares the online performance with the optimal static controller/policy <cit.>. The benchmark in static regret is weaker than that in dynamic regret because the optimal controller may be far from being static, and it has been shown in literature that o(N) static regret can be achieved even without predictions (i.e., W=0). Thus, we will focus on the dynamic regret analysis and study how predictions can improve the dynamic regret. 




[Linear quadratic (LQ) tracking]

	Consider a discrete-time  tracking problem for a system
	x_t+1=  A   x_t +   B   u_t.
	The goal is to minimize the quadratic tracking loss of a trajectory {θ_t}_t=0^N 

    J(  𝐱,   𝐮)=1/2∑_t=0^N-1[ (  x_t-  θ_t)^⊤   Q_t (  x_t-  θ_t)+    u_t^⊤   R_t u_t] + 1/2 (  x_N-  θ_N)^⊤   Q_N(  x_N-  θ_N).
 
	In practice, it is usually difficult to know the complete trajectory {θ_t}_t=0^N a priori, what are revealed are usually the next few steps, making it an online control problem with predictions.
	
	
	



	
	




























Assumptions and  useful concepts.
Firstly, we assume controllability, which is standard in control theory and roughly means that the system can be steered to any state by proper control inputs <cit.>.


	The LTI system x_t+1=Ax_t+Bu_t is controllable.

It is well-known that any controllable LTI system can be linearly transformed to a canonical form <cit.> and the linear transformation can be computed efficiently a priori using A and B, which can further be used to reformulate the cost functions f_t, g_t. Thus, without loss of generality, this paper only considers LTI systems in the canonical form, defined as follows.

A system x_t+1=Ax_t+Bu_t is said to be in the canonical form if 

    A=[
    	[ 0 1   0  ; ⋮ ⋱ ⋱  ;     0 1  ;   * ⋯ * * * … *     ⋯ *;         0 1 0        ;         ⋮ ⋱ ⋱      ;             0 1        ;   * ⋯ * * * ⋯ * ⋯ * ⋯   *; ⋯             ⋯;         0 1 ⋯ 0;         ⋮ ⋱ ⋱;           0 1; * ⋯ * * * ⋯ * ⋯ * * ⋯ * ]],    B= [
    	[  0  0  …;  ⋮  ⋮  ⋮;  0   ;  1  0  ⋯;  0  0   ;  ⋮  ⋮  ⋯;         ;  0  1  ⋯;  ⋯     ⋯;  0  ⋯  ⋯;  ⋮  ⋱;  0  0 ⋯1 ]],

	where each * represents a (possibly) nonzero entry, and the rows of B with 1 are the same rows of A with * and the indices of these rows are denoted as {k_1, …, k_m}. 
	Moreover, let p_i = k_i -k_i-1 for  1 ≤ i ≤ m, where k_0=0. The  of a canonical-form (A,B) is defined as 
    p = max{p_1, …, p_m }.




 












Next, we introduce assumptions on the cost functions and their optimal solutions. 



Assume f_t is μ_f strongly convex and l_f Lipschitz smooth for  0≤ t ≤ N, and g_t is  convex and l_g Lipschitz smooth for 0 ≤ t ≤ N-1 for some μ_f,   l_f, l_g>0. 



	

	
	
	
	
	Assume the minimizers to f_t, g_t, denoted as θ_t = _x f_t(x), ξ_t =  _u g_t(u), are uniformly bounded,  i.e. there exist θ̅, ξ̅ such that θ_t ≤θ̅, ξ_t ≤ξ̅,  ∀  t.
	


	


	

These assumptions are commonly adopted in convex analysis. The uniform bounds  rule out extreme cases. Notice  that the LQ tracking problem in Example 1 satisfies  Assumption <ref> and <ref> if Q_t, R_t are positive definite with uniform bounds on the eigenvalues and if θ_t are uniformly bounded for all t. 






















§ ONLINE CONTROL ALGORITHMS: RECEDING HORIZON GRADIENT-BASED CONTROL






This section introduces our online control algorithms, receding horizon gradient-based control (RHGC). The design is by first converting the online control problem to an equivalent online optimization problem with finite temporal-coupling costs, then designing gradient-based online optimization algorithms by utilizing this finite temporal-coupling property. 

 










 §.§ Problem transformation



Firstly, we notice that the offline optimal control problem (<ref>) can be viewed as an optimization with equality constraints over x and u. The individual stage cost f_t(x_t)+g_t(u_t) only depends on the current x_t and u_t but the equality constraints couple x_t, u_t with x_t+1 for each t. In the following, we will rewrite (<ref>) in an equivalent form of an unconstrained optimization problem on some entries of x_t for all t, but the new stage cost at each time t will depend on these new entries across a few nearby time steps. We will harness this structure to design our online algorithm.  


In particular, the entries of x_t adopted in the reformulation are: x_t^k_1,…, x_t^k_m, where ={k_1, …, k_m} is defined in Definition <ref>. For ease of notation, we define 

    z_t (x_t^k_1,…, x_t^k_m)^⊤,    t≥ 0
 
and write z_t^j=x_t^k_j where j=1,…,m. Let z (z^⊤_1,…, …, z^⊤_N)^⊤. By the  canonical-form equality constraint x_t=Ax_t-1+Bu_t-1, we have x_t^i=x_t-1^i+1 for i∉, so x_t can be represented by z_t-p+1,…, z_t in the following way:

    x_t=(z_t-p_1+1^1, …, z_t^1_p_1, z_t-p_2+1^2, …, z_t^2_p_2, …, z_t-p_m+1^m, …, z_t^m_p_m)^⊤,    t≥0,

where z_t for t≤ 0 is determined by x_0 in a way to let (<ref>) hold for t=0. For  ease of  exposition and without loss of generality,  we consider x_0=0 in this paper; then we have z_t=0 for t≤ 0.
Similarly, u_t can be determined by z_t-p+1, …, z_t, z_t+1 by 

    u_t = z_t+1- A(,:) x_t= z_t+1-A(,:)(z_t-p_1+1^1, …, z_t^1,  …, z_t-p_m+1^m, …, z_t^m)^⊤,    t≥ 0
 
where A(,:) consists of k_1, …, k_m rows of A. 





It is straightforward to verify that equations (<ref>, <ref>, <ref>) describe a bijective transformation between {( x,  u)| x_t+1=Ax_t +Bu_t}  and z ∈^mN, since the LTI constraint x_t+1=Ax_t +Bu_t is naturally embedded in the relation (<ref>, <ref>). Therefore, based on the transformation, an optimization problem with respect to z ∈^m N can be designed to be equivalent with (<ref>). Notice that the resulting optimization problem  has no constraint on z. Moreover, the cost functions on z can be obtained by substituting  (<ref>, <ref>) into f_t(x_t) and g_t(u_t), i.e. f̃_t(z_t-p+1, …, z_t) f_t(x_t) and g̃_t(z_t-p+1, …, z_t, z_t+1) g_t(u_t). 
Correspondingly, the objective function of the equivalent optimization with respect to z is

    C(𝐳)∑_t=0^Nf̃_t(z_t-p+1,…, z_t)+∑_t=0^N-1g̃_t(z_t-p+1,…,z_t+1)

 C( z) has many nice properties, some of which are formally stated  below. 
 
The function C( z) has the following properties:



-1pt
    
  i) C( z) is μ_c=μ_f strongly convex and l_c smooth for  l_c= pl_f + (p+1)l_g [I_m,-A(,:)]^2.
    
  ii) For any ( x, u) s.t. x_t+1= Ax_t+B u_t, C( z)=J( x, u) where z is defined in (<ref>). Conversely,   ∀  z, the ( x,  u) determined by (<ref>,<ref>) satisfies   x_t+1= A x_t+ B u_t and J( x, u)=C( z);
    
  iii) Each stage cost f̃_t+g̃_t in (<ref>) only depends on z_t-p+1,…,z_t+1. 



Property ii) implies that any online algorithm for deciding z can be translated to an online algorithm for x and u by (<ref>, <ref>) with the same costs.

Property iii) highlights one nice property, finite temporal-coupling, of C( z), which serves as a foundation for our online algorithm design.  







 For illustration,  
consider the following dynamical system with n=2, m=1:

    [
    	[ x_t+1^1; x_t+1^2 ]]= [
    	[   0   1; a_1 a_2 ]] [
    	[ x_t^1; x_t^2 ]]+ 
    	[
    	[ 0; 1 ]]u_t

	Here, k_1=2, ℐ={2}, A(ℐ,:)=(a_1, a_2), and z_t=x_t^2.  By (<ref>),
 x_t^1= x_t-1^2 and x_t=(z_t-1,z_t)^⊤. Similarly, u_t=x^2_t+1- A(ℐ,:) x_t=z_t+1- A(ℐ,:)(z_t-1,z_t)^⊤. Hence,  f̃_t(z_t-1, z_t)= f_t(x_t)=f_t((z_t-1,z_t)^⊤), g̃_t(z_t-1,z_t, z_t+1)= g_t(u_t)=g_t(z_t+1- A(ℐ,:) (z_t-1,z_t)^⊤ ).  

 This paper considers a  reparameterization method with respect to states x via the canonical form, and it might be interesting to compare it with the  more direct reparameterization with respect to control inputs u. The control-based reparameterization has been discussed in literature <cit.>. It has been observed in <cit.> that when A is not stable, the condition number of the cost function derived from the control-based reparameterization goes to infinity as W→ +∞, which may result in computation issues when W is large. However, the state-based reparameterization considered in this paper can guarantee bounded condition number for all W even for unstable A, as shown in Lemma <ref>. This is one major advantage of the state-based reparameterization method considered in this paper.







 §.§ Online algorithm design: RHGC

This section introduces our RHGC  based on the reformulation (<ref>) and inspired by  <cit.>. As mentioned earlier, any online algorithm for z_t can be translated to  an online algorithm for x_t,u_t. Hence, we will focus on designing an online algorithm for z_t in the following. By the finite temporal-coupling property of C( z), the partial gradient of the total cost C( z) only depends on the finite neighboring stage costs {f̃_τ, g̃_τ}_τ=t^t+p-1 and finite neighboring stage variables (z_t-p,…,z_t+p)=:z_t-p:t+p. 

    ∂ C/∂ z_t( z)= ∑_τ=t^t+p-1∂f̃_τ/∂ z_t(z_τ-p+1,…, z_τ) + ∑_τ=t-1^t+p-1∂g̃_τ/∂ z_t(z_τ-p+1,…, z_τ+1)

Without causing any confusion, we use ∂ C/∂ z_t(z_t-p:t+p) to denote ∂ C/∂ z_t( z) for highlighting the local dependence. Thanks to the local dependence, despite  the fact that not all the future costs are available, it is still possible to compute the partial gradient of the total cost by using only a finite lookahead window of the cost functions. This observation motivates the design of our receding horizon gradient-based control (RHGC) methods, which are the online implementation of gradient methods, such as  vanilla gradient descent, Nesterov's accelerated gradient, triple momentum, etc., <cit.>. 
 
 









Firstly, we illustrate the main idea of RHGC by receding horizon gradient descent (RHGD) based on vanilla gradient descent. In RHGD (Algorithm <ref>), index j refers to the iteration number of the corresponding gradient update of C( z). There are two major steps to decide z_t. Step 1 is initializing the decision variables z(0). Here, we do not restrict the initialization algorithm φ and allow any oracle/online algorithm without using lookahead information, i.e.  z_t+W(0) is selected based only on the information up to t+W-1: z_t+W(0)=φ({f̃_s, g̃_s }_s=0^t+W-1).
One example of φ will be provided in Section 4. Step 2 is using the W-lookahead costs to conduct gradient updates. Notice that the gradient update from z_τ(j-1) to z_τ(j) is implemented in a backward order of τ, i.e. from τ=t+W to τ=t. Moreover, since the partial gradient  ∂ C/∂ z_t requires the local decision variables z_t-p:t+p, given W-lookahead information,  RHGD can only conduct K=W-1/p iterations of gradient descent for the total cost  C( z). For more  discussion,  we  refer the reader to <cit.> for the p=1 case. 



In addition to RHGD, RHGC can also incorporate accelerated gradient methods in the same way, such as Nesterov's accelerated gradient and triple momentum. For the space limit, we only formally present receding horizon triple momentum (RHTM) in Algorithm <ref> based on triple momentum <cit.>. 
RHTM also consists of two major steps when determining z_t: initialization and gradient updates based on the lookahead window. The two major differences from RHGD are that the decision variables in RHTM include not only z(j) but also auxiliary variables ω(j) and y(j), which are adopted in triple momentum to accelerate the convergence, and  that the gradient update  is by triple momentum instead of gradient descent. Nevertheless, RHTM can also conduct K=W-1/p iterations of triple momentum for  C( z) since the triple momentum update requires the same neighboring cost functions.











Though it appears that RHTM does not fully exploit the lookahead information since only a few gradient updates are used, in Section <ref>, we show that RHTM  achieves near-optimal performance with respect to W, which means that RHTM successfully extracts and utilizes the prediction information. 



Finally, we briefly introduce MPC<cit.> and suboptimal MPC<cit.>, and compare them with our algorithms. MPC tries to solve a W-stage optimization at each  t and implements the first control input. Suboptimal MPC, as a variant of MPC aiming at reducing computation, conducts an optimization method only for a few iterations without solving the optimization completely. Our algorithm's computation time is similar to that of suboptimal MPC with a few gradient iterations. 
However, the major difference between our algorithm and suboptimal MPC is that suboptimal MPC conducts gradient updates for a truncated W-stage optimal control problem based on W-lookahead information, while our algorithm is able to conduct  gradient updates for the complete N-stage optimal control problem based on the same W-lookahead information by utilizing the reformulation  (<ref>, <ref>, <ref>, <ref>). 








 







































































	

	
	

	


























 

























	

















































































§ REGRET UPPER BOUNDS


Because our RHTM (RHGD) is designed to exactly implement the triple momentum (gradient descent) of C( z) for K iterations, it is straightforward to have the following regret guarantees  that connect  the regrets of RHTM and RHGD with the regret of the initialization oracle φ,





 
	Consider W≥ 1 and stepsizes γ_g=1/l_c, γ_c = 1+ϕ/l_c, γ_ω = ϕ^2/2-ϕ, γ_y = ϕ^2/(1+ϕ)(2-ϕ), γ_z = ϕ^2/1-ϕ^2,  ϕ= 1-1/√(ζ), 
	and let ζ=l_c/μ_c denote C( z)'s condition number. For any oracle φ, 
	
    (RHGD) ≤ζ(ζ-1/ζ)^K(φ),   (RHTM) ≤ζ^2 (√(ζ)-1/√(ζ))^2K(φ)

	where K= W-1/p, (φ) is the regret of  the initial controller: u_t(0)=z_t+1(0)-A(, :)x_t(0).








 Theorem <ref> suggests that  for any online algorithm φ without predictions, RHGD and RHTM  can use predictions to lower the regret by a factor of ζ (ζ-1/ζ)^K and ζ^2(√(ζ)-1/√(ζ))^2K respectively via additional K=W-1/p gradient updates.
Moreover, the factors decay exponentially with K=W-1/p, and K almost linearly increases with W. This indicates that  RHGD and RHTM   improve the performance exponentially fast with an increase in the prediction window W for any initialization method. In addition, K=W-1/p decreases with  p, implying that the regrets increase with the controllability index p (Definition <ref>). This is intuitive because p roughly indicates how fast the controller can influence the system state effectively: the larger the p is, the longer it takes. 
To see this, consider Example <ref>. Since u_t-1 does not directly affect x_t^1, it takes at least p=2 steps to change x_t^1 to a desirable value. Finally, RHTM's regret decays faster than RHGD's, which is intuitive because triple momentum converges faster than gradient descent. Thus, we will focus on  RHTM in the following.






An initialization method: follow the optimal steady state (FOSS).
To complete the regret analysis for RHTM, we provide a simple initialization method, FOSS, and its dynamic regret bound.


As
mentioned before, any online control algorithm without predictions, e.g.  <cit.>, can be applied as an initialization oracle φ. However, most literature  study  static regrets  rather than dynamic regrets.




The optimal steady state for stage cost f(x)+g(u) refers to (x^e,u^e)_x=Ax+Bu(f(x)+g(u)).

 Follow the optimal steady state algorithm (FOSS) first solves the optimal steady state (x_t^e, u_t^e) for cost  f_t(x)+g_t(u), then determines z_t+1 by x_t^e, i.e. 
 z_t+1=(x_t^e,k_1, …, x_t^e,k_m)^⊤ at each t+1.

FOSS is motivated by the fact that the optimal steady state cost is the optimal infinite-horizon average cost for LTI systems with time-invariant cost functions <cit.>, so FOSS should yield acceptable performance at least for slowly changing cost functions. Nevertheless, we admit that  FOSS is proposed mainly for analytical purposes and other online algorithms may outperform FOSS in various perspectives.
 The following is a regret bound for FOSS, relying on the solution to  Bellman equations.

Consider optimal control problem: minlim_N→ +∞1/N∑_t=0^N-1 (f(x_t)+g(u_t)) where x_t+1=Ax_t + Bu_t.
    Let λ^e be the optimal steady state cost f(x^e)+g(u^e), which is also the optimal infinite-horizon average cost <cit.>. The Bellman equations for the    problem is h^e(x)+λ^e = min_u (f(x)+g(u)+h^e(Ax+Bu)).



    The solution to the Bellman equations, denoted by h^e(x), is sometimes called as a bias function <cit.>. To ensure the uniqueness of the solution,  some extra conditions, e.g. h^e(0)=0, are usually  imposed.
   













Let (x_t^e,u_t^e) and h_t^e(x) denote the optimal steady state and the bias function with respect to cost f_t(x)+g_t(u) respectively for 0≤ t ≤ N-1. Suppose h_t^e(x) exists for 0≤ t ≤ N-1,[h_t^e may not be unique, so extra conditions  can be imposed on h_t^e for more interesting    regret bounds.] then the regret of FOSS can be bounded by

    (FOSS)= O(∑_t=0^N( x_t-1^e- x_t^e  +h^e_t-1(x_t^*)-h^e_t(x_t^*))),

	where
	  {x_t^*}_t=0^N denotes the optimal state trajectory for (<ref>), x_-1^e=x_0^*=x_0=0, h^e_-1(x)=0, h^e_N(x)=f_N(x), x_N^e=θ_N.
	  Consequently, by Theorem <ref>, the regret bound of RHTM with initialization FOSS is
(RHTM)= O(ζ^2(√(ζ)-1/√(ζ))^2K∑_t=0^N( x_t-1^e- x_t^e  +h^e_t-1(x_t^*)-h^e_t(x_t^*))).


Theorem <ref> bounds the regret by the variation of the optimal steady states x_t^e and the bias functions h_t^e. 
If f_t and g_t do not change, x_t^e and  h_t^e do not change, yielding a small O(1) regret, i.e. O(x_0^e+h_0^e(x_0)),  matching our intuition. Though Theorem <ref> requires  h_t^e exists, the existence is guaranteed for many control problems, e.g. LQ tracking and control problems with turnpike properties <cit.>. 


















§ LINEAR QUADRATIC TRACKING: REGRET UPPER BOUNDS AND A FUNDAMENTAL LIMIT


To provide more intuitive meaning for our regret analysis in Theorem <ref> and Theorem <ref>, we apply RHTM to the LQ tracking problem in Example <ref>. Results for the time varying Q_t, R_t, θ_t are provided in  Appendix E; whereas here we focus on a special case which gives clean expressions for regret bounds: both an upper bound for RHTM with initialization FOSS and a lower bound for any online algorithm. Further, we show that the lower bound and the upper bound almost match each other, implying that our online algorithm RHTM uses the predictions in a nearly optimal way even though it only conducts a few gradient updates at each time step . 

The special case of LQ tracking problems is in the following form,
	
    1/2∑_t=0^N-1[ (  x_t-  θ_t)^⊤   Q (  x_t-  θ_t)+    u_t^⊤   R u_t] + 1/2   x_N^⊤   P^e  x_N,

	where Q>0, R>0, and P^e is the solution to the algebraic Riccati equation with respect to Q, R <cit.>. Basically, in this special case, Q_t=Q, R_t=R for 0≤ t ≤ N-1, Q_N=P^e, θ_N=0, and only θ_t changes for t=0, 1, …, N-1. The  LQ  tracking problem (<ref>) aims to follow a time-varying trajectory {θ_t} with constant weights on the tracking cost and the control cost. 







Regret upper bound. Firstly, based on Theorem <ref> and Theorem <ref>, we have the following bound.

	
	Under the stepsizes in Theorem <ref>,  RHTM with FOSS as the initialization rule satisfies
	
    (RHTM) = O(ζ^2 (√(ζ)-1/√(ζ))^2K∑_t=0^Nθ_t-θ_t-1)

	where  K= (W-1)/p,   ζ is the condition number of the corresponding C( z),  θ_-1=0.	

This corollary shows that the regret can be bounded by the total variation of θ_t for constant Q, R. 

Fundamental limit. For any online algorithm, we have the following lower bound.


Consider 1 ≤ W≤ N/3, any condition number ζ >1, any variation budget 4 θ̅≤ L_N ≤  (2N+1) θ̅, and any controllability index p≥ 1. For any online algorithm , there exists an LQ tracking problem in form (<ref>) where i) the canonical-form system (A, B) has controllability index p, ii) the sequence {θ_t} satisfies the variation budget ∑_t=0^N θ_t-θ_t-1≤ L_N, and iii) the corresponding C( z) has condition number ζ, such that the following  lower bound holds
	
    (𝒜) = Ω(  (√(ζ)-1/√(ζ)+1)^2 KL_N )= Ω(  (√(ζ)-1/√(ζ)+1)^2 K∑_t=0^N θ_t-θ_t-1)

	where   K = (W-1)/p and  θ_-1=0.

Firstly, the lower bound in Theorem <ref> almost matches the upper bound in Corollary <ref>, especially when ζ is large, demonstrating that RHTM utilizes the predictions  in a near-optimal way. The major conditions in Theorem <ref> require that the prediction window is short compared with the horizon: W≤ N/3, and the variation of the cost functions should not be too small: L_N≥ 4θ̅, otherwise the online control problem is too easy and the regret can be very small. Moreover, the small gap between the regret bounds is conjectured to be nontrivial, because this gap coincides with the long lasting  gap in the convergence rate of the first-order algorithms for strongly convex and smooth optimization. In particular, the lower bound in Theorem <ref> matches the fundamental convergence limit  in <cit.>, and the upper bound is by triple momentum's convergence rate, which is the best one to our knowledge. 









§ NUMERICAL EXPERIMENTS
















LQ tracking problem in Example <ref>. The system considered here has n=2, m=1, and p=2. More details  of the experiment settings are provided in Appendix H. 
 We compare RHGC with a suboptimal MPC algorithm, fast gradient MPC (subMPC) <cit.>. Roughly speaking, subMPC solves the W-stage truncated optimal control from t to t+W-1 by Nesterov's accelerated gradient <cit.>, and one iteration of Nesterov's accelerated gradient requires 2W gradient evaluations of stage cost function since W stages are considered and each stage has two costs f_t and g_t.  This implies that, in terms of the number of gradient evaluations, subMPC with one iteration corresponds to our RHTM because RHTM also requires roughly 2W gradient evaluations per stage. Therefore, Figure <ref> compares our RHGC algorithms with subMPC with one iteration. Figure <ref> also plots subMPC with 3 and 5 iterations for more insights. Besides, Figure <ref> plots not only RHGD and RHTM, but also RHAG, which is based on Nesterov's accelerated gradient.  
Figure <ref> shows that all our  algorithms achieve exponential decaying regrets with respect to W, and the regrets are piecewise constant, matching Theorem <ref>.  Further, it is observed that RHTM and RHAG perform better than RHGD, which is intuitive because triple momentum and Nesterov's accelerated gradient are accelerated versions of gradient descent. In addition, our algorithms are much better than   subMPC with 1 iteration, implying that our algorithms utilize the lookahead information more efficiently given  similar computational time. Finally,
 subMPC achieves better performance by increasing the iteration number but the improvement saturates as W increases, in contrast to the steady improvement of RHGC.  






Path tracking for a two-wheel mobile robot. Though we presented our online algorithms on an LTI system, our RHGC methods are applicable to some nonlinear systems as well. Here we consider a two-wheel mobile robot with nonlinear kinematic dynamics  ẋ = vcosδ, ẏ = v sinδ, δ̇= w



where (x,y) is the robot location, v and  w are the tangential and angular velocities respectively, δ denotes the tangent angle between v and the x axis <cit.>. The control is directly on the v and w, e.g., via the pulse-width modulation (PWM) of the motor <cit.>.
Given a reference path (x^r_t, y^r_t), the objective is to balance the tracking performance and the control cost, i.e., min  ∑_t = 0^N [ c_t·((x_t - x^r_t)^2 + (y_t - y^r_t)^2 )+c^v_t· (v_t)^2 + c^w_t· (w_t)^2].





We discretize the dynamics with time interval Δ t=0.025s; then  follow similar ideas in this paper to reformulate the optimal path tracking problem  to an unconstrained optimization with respect to (x_t,y_t) and apply RHGC. See Appendix H for details. Figure <ref> plots the tracking results with window W=40 and W=80 corresponding to lookahead time 1s and 2s. A video showing the dynamic processes with different W is provided at <https://youtu.be/fal56LTBD1s>. It is observed that the robot follows the reference trajectory well especially when the path is smooth but deviates a little more when the path has sharp turns, and a longer lookahead window leads to better tracking performance.  These results confirm that our RHGC works effectively on nonlinear systems.




























§ CONCLUSION

This paper studies the role of predictions on dynamic regrets of online control problems with linear dynamics. We design RHGC algorithms and provide  regret upper bounds of two specific algorithms: RHGD and RHTM. We also provide a fundamental limit and show the fundamental limit almost matches RHTM's upper bound. This paper leads to many interesting future directions, some of which are briefly discussed below. The first direction is to study more realistic prediction models which considers random prediction noises, e.g. <cit.>.  The second direction is to consider unknown systems with process noises, possibly by applying learning-based control tools <cit.>. Further, more studies could be conducted on general control problems including nonlinear control and control with input and state constraints. Besides, it is interesting to consider other performance metrics, such as competitive ratio, since the dynamic regret is non-vanishing. Finally, other future directions include closing the gap of the regret bounds and more discussion on the effect of the canonical-form transformation on the condition number.














§ ACKNOWLEDGEMENT

This work was supported by NSF Career 1553407, ARPA-E NODES, AFOSR YIP and ONR YIP programs. 


unsrt




§ APPENDICES


In Appendix <ref>, we will discuss the canonical-form transformation. In Appendix <ref>, we will briefly introduce the  triple momentum algorithm proposed in <cit.> and provide the proof of Theorem <ref>. In Appendix <ref>, we will provide the proof of Lemma <ref>. In Appendix <ref>, we will present the proof of Theorem <ref>.  Appendix <ref> provides a proof of Corollary <ref> and regret analysis for more general linear quadratic tracking problems. Appendix <ref> provides a proof of Theorem <ref>. In Appendix <ref>, we will provide the proofs of the technical lemmas used in Appendix <ref>. In Appendix <ref>, we will provide  a more detailed description of our numerical experiments.



§ CANONICAL FORM


In this section, we introduce the linear transformation from a general LTI system to a canonical-form LTI system, and then discuss how to convert a general online optimal control problem to an online optimal control problem with a canonical-form system.


Firstly, consider a general LTI system: x_t+1=Ax_t +Bu_t and two invertible matrices S_x∈^n, S_u∈ R^m. Under the linear transformation on state and control: x̂_t= S_x x_t,  û_t = S_u u_t, the equivalent LTI system with respect to the new state x̂_t and the new control û_t is

    x̂_t+1= S_x A S_x^-1x̂_t + S_x B S_u^-1û_t


By   <cit.>, for any controllable (A, B), there exist S_x, S_u such that Â = S_x A S_x^-1 and B̂= S_x B S_u^-1 are in the canonical form defined in Definition <ref>. The computation method of S_x, S_u is also provided in <cit.>. 


In an online optimal control problem, since A,B are known as priors, S_x,  S_u can be computed offline. When stage cost functions f_t(x_t), g_t(u_t) are received online, the new cost functions f̂_t(x̂_t), ĝ_t(û_t) for the canonical-form system can be computed online by applying S_x,S_u:

    f̂_t(x̂_t)= f_t(x_t)=f_t(S_x^-1x̂_t),  ĝ_t(û_t)= g_t(u_t)=g_t(S_u^-1û_t)

Moreover, it is straightforward to verify that f̂_t(x̂_t) and ĝ_t(û_t) still satisfy Assumption <ref> and <ref>, just with perhaps different parameters. For example, f̂_t(x̂_t) is μ_f/S_x^2 strongly convex and l_f S_x^-1^2 smooth and ĝ_t(û_t) is  convex  and l_g S_u^-1^2 smooth.   Therefore, it is without loss of generality to only consider online optimal control with canonical-form systems.








§ TRIPLE MOMENTUM AND A PROOF OF THEOREM <REF>




Triple Momentum (TM) is an accelerated version of gradient descent proposed in <cit.>. When optimizing an unconstrained optimization min_ zC( z), at each iteration j≥ 0, TM conducts

    ω(j+1)= (1+δ_ω) ω(j)-δ_ωω_(j-1) -δ_c ∇ C( y(j))
        y(j+1) = (1+δ_y) ω(j+1) -δ_yω(j)
        z(j+1)=(1+δ_z)ω(j+1)-δ_z ω(j)

where ω(j),  y(j) are auxiliary variables to accelerate the convergence, z(j) is the decision variable, ω(0)=ω(-1)= z(0)= y(0) are given initial values. 



Suppose z=(z_1^⊤ , …, z_N^⊤ )^⊤. Zooming in to each coordinate z_t, the update of z_t(j) by TM is provided below

    ω_t(j+1)= (1+δ_ω) ω_t(j)-δ_ωω_t(j-1) -δ_c  ∂ C/∂ y_t( y(j))
        y_t(j+1) = (1+δ_y) ω_t(j+1) -δ_y ω_t(j)
        z_t(j+1)=(1+δ_z)ω_t(j+1)-δ_z ω_t(j)

By Section <ref>,  ∂ C/∂ y_t( y(j)) only depends on stage cost functions and  stage variables  across a finite neighboring stages, 
allowing the online implementation in Algorithm <ref> based on the finite-lookahead window.



TM enjoys a faster convergence rate than the gradient descent for μ_c strongly convex and l_c smooth functions under proper step sizes. In particular, when γ_c = 1+ϕ/l_c,  γ_w = ϕ^2/2-ϕ,  γ_y = ϕ^2/(1+ϕ)(2-ϕ),  γ_z = ϕ^2/1-ϕ^2, and ϕ= 1-1/√(ζ), ζ= l_c/μ_c, by <cit.>, the convergence of TM satisfies:

    C( z(j))-C( z^*) ≤ (√(ζ)-1/√(ζ))^2jl_c ζ/2 z(0)- z^* ^2 ≤ζ^2 (√(ζ)-1/√(ζ))^2j(C(𝐳(0))-C(𝐳^*))

In the following, we will apply this result to prove Theorem <ref>.










 §.§ Proof of Theorem <ref>



By comparing TM with RHTM, it can be verified that z_t+1(K) computed by RHTM is the same as  z_t+1(K) computed by the triple momentum after  K iterations. Moreover,
	according to  the equivalence between  the optimization  min_ z C( z) and the optimal control J( x,  u) in Lemma <ref>,  
    J(RHTM)=C( z(K)),  J(φ)=C( z(0)),  J^*=C( z^*)

	Finally, by utiltizing (<ref>), the bound on Regret(RHTM) is straightforward.


The regret of RHGD can be proved in the same way.









§ PROOF OF LEMMA <REF>

Property ii) and iii) can be directly verified by definition. Thus, it suffices to prove i): the  strong convexity and smoothness of C( z).

Notice that x_t,  u_t are linear with respect to z by (<ref>) (<ref>). For ease of reference, we define matrix M^x_t, M^u_t to represent the relation between x_t,u_t and z, i.e,  x_t = M^x_t z and u_t = M^u_t z. Similarly, we write f̃_t(z_t-p+1, …, z_t) and g̃_t(z_t-p+1, …, z_t+1) in terms of z for simplicity of notation:

    f̃_t(z_t-p+1, …, z_t)   = f̃_t( z) = f_t(M^x_t z)
    g̃_t(z_t-p+1, …, z_t+1)   = g̃_t( z) = g_t(M^u_t z)


A direct consequence of the linear relations is that f̃_t( z) and g̃_t( z) are convex with respect to z because f_t(x_t), g_t(u_t) are convex and the linear transformation preserves convexity. 

In the following, we will focus on the proof of strong convexity and smoothness. For simplicity, in the following, 	 we  only consider cost function f_t, g_t with minimum values  zero: f_t(θ_t)=0, and g_t(ξ_t)=0 for all t. This is  without loss of generality because by strong convexity and smoothness, f_t, g_t have minimum values, and by subtracting the minimum value, we can let f_t, g_t have minimum value 0.



  
Strong convexity. Since g̃_t is convex, we only need to prove that ∑_t f̃_t( z) is strongly convex then the sum C( z) is strongly convex because the sum of convex functions and a strongly convex function is strongly convex.

In particular, by the strong convexity of f_t(x_t), we have the following result: for any z,  z'∈^Nm and x_t=M^x_t z, x_t'=M^x_t z':
	
    f̃_t( z')   - f̃_t( z) - ⟨∇f̃_t( z),  z' - z ⟩ - μ_f/2 z'_t - z_t ^2
        = f̃_t( z')- f̃_t( z) - ⟨ (M^x_t)^⊤∇ f_t(x_t),  z' - z ⟩ - μ_f/2 z'_t - z_t ^2
        =  f̃_t( z')- f̃_t( z) - ⟨∇ f_t(x_t), M^x_t( z' - z) ⟩ - μ_f/2z'_t - z_t ^2
       = f̃_t( z')- f̃_t( z) - ⟨∇ f_t(x_t),  x'_t -x_t ⟩ - μ_f/2 z'_t - z_t ^2   
       ≥ f_t(x'_t)-f_t(x_t )- ⟨∇ f_t(x_t),  x'_t -x_t⟩ - μ_f/2 x'_t - x_t ^2  ≥ 0

where the first equality is by the chain rule, the second equality is by the definition of inner product, the third equality is by the definition of x_t, x_t', the first inequality is by f̃_t(z)=f_t(x) and z_t=(x_t^k_1,…, x_t^k_m)^⊤, and the last inequality is because
f_t(x_t) is μ_f strongly convex.

Summing over t on both sides of the inequality results in the strong convexity of ∑_t f̃_t( z):

    ∑_t=1^N[f̃_t( z')- f̃_t( z) - ⟨∇f̃_t( z),  z' - z ⟩ - μ_f/2 z'_t - z_t ^2 ]  
    
    	=   ∑_t=1^Nf̃_t( z')- ∑_t=1^Nf̃_t( z) - ⟨∇∑_t=1^Nf̃_t( z),  z' - z ⟩ - μ_f/2 z' -  z ^2 ≥ 0

	Consequently, C( z) is  strongly convex with parameter at least μ_f  by the convexity of g̃_t.





  
Smoothness. 
We will prove the smoothness by considering  f̃_t( z) and g̃_t( z) respectively. 
	
	Firstly, let's consider f̃_t( z). Similar to the proof for strong convexity, we use the smoothness of f_t(x_t). For any z,  z', and x_t=M^x_t z, x_t'=M^x_t z',  we can show that
	
    f̃_t( z')    = f_t( x'_t) ≤ f_t(x_t)+ ⟨∇ f_t(x_t),  x'_t-x_t⟩ + l_f/2 x'_t-x_t^2
       ≤f̃_t( z) + ⟨∇f̃_t( z),  z'- z ⟩ + l_f/2( z'_t-p+1-z_t-p+1^2 + … +  z'_t-z_t^2 )

where the second inequality is by x_t=M^x_t z and the chain rule and (<ref>). 
	
	Secondly, we consider g̃_t(z) in a similar way. For any z,  z', and u_t=M^u_t z, u_t'=M^u_t z', we have
	
    g̃_t( z')	   = g_t( u'_t)≤ g_t(u_t)+  ⟨∇ g_t(u_t),  u'_t-u_t⟩ + l_g/2u'_t- u_t^2 
        =g̃_t( z)+  ⟨  (M^u_t)^⊤∇ g_t(u_t),  z'- z⟩ + l_g/2u'_t- u_t^2  
        = g̃_t( z)+  ⟨∇g̃_t( z),  z- z⟩ + l_g/2u'_t- u_t^2

		Since u_t= z_t+1-A(, :) x_t=[I_m,-A(, :)](z_t+1^⊤ ,x_t^⊤ )^⊤, we have that 
		
    l_g/2u'_t- u_t^2  	   ≤l_g/2[I_m,-A(, :)][ ( (z_t+1')^⊤ , (x'_t)^⊤ )^⊤  - ( z_t+1^⊤ ,  x^⊤ _t)^⊤] ^2 
       ≤l_g/2[I_m,-A(, :)]^2 (z_t+1- z'_t+1^2+ x_t -  x'_t^2)
       ≤l_g/2[I_m,-A(, :)]^2 (z_t+1-z'_t+1^2+ … + z_t-p+1-z'_t-p+1^2)

	Finally, by summing f̃_t( z'), g̃_t( z')'s inequalities above over all t, we have
		
    C( z') ≤ C( z) + ⟨∇ C( z),  z'- z⟩ + (p l_f + (p+1)l_g [I_m,-A(, :)]^2)/2  z'- z^2

 Thus, we have proved the smoothness of C( z).
	
	
	
	



























§ PROOF OF THEOREM <REF>

Remember that Regret(FOSS) =  J(FOSS)-J^*. To  bound the regret, we let the sum of the optimal steady state costs, ∑_t=0^N-1λ^e_t, be a middle ground and bound J(FOSS)-∑_t=0^N-1λ^e_t and 
∑_t=0^N-1λ^e_t-J^* in Lemma <ref> and Lemma <ref> respectively. Then, the regret bound can be obtained by combining the two bounds.

 Let  x_t(0) denote the state determined by FOSS. 

    J(FOSS)-∑_t=0^N-1λ^e_t ≤ c_1 ∑_t=0^N-1x_t^e -x_t-1^e + f_N(x_N(0))=O(  ∑_t=0^Nx_t^e -x_t-1^e)

   where we define x_N^eδ_N, x_-1^e x_0=0 for simplicity of  notation, c_1 is a constant that does not depend on N, W and big O hides a constant that does not depend on N, W.




Let h_t^e(x) denote a solution to the  Bellman equations under cost f_t(x)+g_t(u). Let { x_t^*} denote the optimal state trajectory to the offline optimal control (<ref>). 

    ∑_t=0^N-1λ^e_t-J^*≤∑_t=1^N(h^e_t-1(x_t^*)-h^e_t(x_t^*))-h_0^e(x_0) = ∑_t=0^N(h^e_t-1(x_t^*)-h^e_t(x_t^*))

   where we define h_N^e(x) f_N(x), h_-1^e(x) 0 and x_0^* x_0 for simplicity of  notation. 


Then, we can complete the proof by applying 
 Lemma <ref> and <ref>:

    J(FOSS)-J^*    =J(FOSS)-∑_t=0^N-1λ^e_t + ∑_t=0^N-1λ^e_t-J^*
       = O(∑_t=0^N( x_t-1^e- x_t^e   + h^e_t-1(x_t^*)-h^e_t(x_t^*)))


In the following, we will prove Lemma <ref> and <ref> respectively. For simplicity, 	 we  only consider cost function f_t, g_t with minimum values zero: f_t(θ_t)=0, and g_t(ξ_t)=0 for all t. There is no loss of generality because by strong convexity and smoothness, f_t, g_t have minimum values, and by subtracting the minimum value, we can let f_t, g_t have minimum value 0.



 §.§ Proof of Lemma <ref>.




Notice that J(FOSS)= ∑_t=0^N-1(f_t(x_t(0))+g_t(u_t(0)) ) + f_N(x_N(0)) and ∑_t=0^N-1λ_t^e =∑_t=0^N-1(f_t(x_t^e)+g_t(u_t^e) ). Thus, it suffices to bound f_t(x_t(0))-f_t(x_t^e) and g_t(u_t(0))-g_t(u_t^e) for 0≤ t ≤ N-1. We will first focus on f_t(x_t(0))-f_t(x_t^e), then bound g_t(u_t(0))-g_t(u_t^e) in the same way.

For 0≤ t≤ N-1, by the convexity of f_t, and the property of L_2 norm,

    f_t(x_t(0))-f_t(x_t^e)   ≤⟨∇ f_t(x_t(0)), x_t(0)-x_t^e⟩≤∇ f_t(x_t(0))  x_t(0)-x_t^e

In the following, we will bound ∇ f_t(x_t(0)) and x_t(0)-x_t^e respectively. 

Firstly, we provide a bound on ∇ f_t(x_t(0)):

    ∇ f_t(x_t(0)) = ∇ f_t(x_t(0)) -∇ f_t(θ_t)≤ l_f  x_t(0)-θ_t≤ l_f(√(n)x̅^e + θ̅)

where the first equality is because θ_t is the global minimizer of f_t, and first inequality is by Lipschitz smoothness, the second inequality is by θ_t ≤θ̅ according to Assumption <ref> and by x_t(0)≤√(n)x̅^e proved in the following lemma. 


	There exist x̅^e and u̅^e that are independent of N, W, such that x_t^e≤x̅^e and u_t^e ≤u̅^e for all 0≤ t≤ N-1. Moreover, x_t(0)≤√(n)x̅^e for 0≤ t ≤ N and  u_t(0)≤√(n)u̅^e for 0≤ t ≤ N-1, where  x_t(0), u_t(0) denote the state and control at t determined by FOSS.

The proof is technical and is deferred to Appendix <ref>. 



Secondly, we provide a bound on x_t(0)-x_t^e. The proof relies on the expressions of the steady state x_t^e and the initialized state x_t(0) of a canonical-form system.


    Consider a canonical-form system: x_t+1=Ax_t+B u_t. 
    
        
  (a) Any steady state (x,u) is in the form of 
        
    x   =(z^1, …, z^1_p_1, z^2, …, z^2_p_2, …, z^m, …, z^m_p_m)^⊤
    
    u    = (z^1, …, z^m)^⊤ -A(, :) x

for some z^1, …, z^m ∈.
Let z=(z^1, …, z^m)^⊤. 
For the optimal steady state with respect to cost f_t+g_t, we denote the corresponding z as z_t^e, and the optimal steady state can be represented as
 x_t^e=(z_t^e,1, …, z_t^e,1, z_t^e,2, …, z_t^e,2, …, z_t^e,m, …, z_t^e,m)^⊤ and u_t^e= z_t^e - A(, :) x_t^e for 0≤ t ≤ N-1. 

  (b) By FOSS initialization, z_t+1(0)=z_t^e, and x_t(0), u_t(0) satisfy

    x_t(0)   =(z_t-p_1^e,1, …, z_t-1^e,1_p_1, z_t-p_2^e,2,…, z_t-1^e,2_p_2, …, z_t-p_m^e,m,…, z_t-1^e,m_p_m),     0 ≤ t ≤ N
    
        u_t(0)     = z_t^e -A(, :) x_t(0)    0 ≤ t ≤ N-1

where z_t^e=0 for t≤ -1.



    
        
  (a) This is by the definition of the canonical form and the definition of the steady state.
        
  (b) By the initialization, z_t(0)=x^e,_t-1= z_t-1^e. By the relation between z_t(0) and x_t(0), u_t(0), we have x_t^(0)= z_t(0)=z_t-1^e, and x_t^-1(0)=z_t-1(0)=z_t-2^e, so on and so forth. This proves the structure of x_t(0). The structure of u_t(0) is  because u_t(0)=z_t+1(0)-A(, :) x_t(0)= z_t^e -A(, :) x_t(0)
    





By Lemma <ref>, we can bound x_t(0)-x_t^e for 0≤ t ≤ N-1 by 

    x_t(0)-x_t^e    ≤√(z_t-1^e-z_t^e^2 + … +  z_t-p^e - z_t^e^2)
       ≤√(x_t-1^e-x_t^e^2 + … +  x_t-p^e - x_t^e^2)
       ≤x_t-1^e-x_t^e + … +  x_t-p^e - x_t^e
       ≤ p(x_t-1^e-x_t^e + … +  x_t-p^e - x_t-p+1^e)



Combining (<ref>) (<ref>) and (<ref>) yields

    ∑_t=0^N-1 f_t(x_t(0))-f_t(x_t^e)   ≤∑_t=0^N-1∇ f_t(x_t(0))  x_t(0)-x_t^e 
       ≤∑_t=0^N-1l_f(√(n)x̅^e+θ̅) p(x_t-1^e-x_t^e + … +  x_t-p^e - x_t-p+1^e) 
       ≤ p^2 l_f(√(n)x̅^e+θ̅) ∑_t=0^N-1x_t-1^e-x_t^e

Notice that the constant term p^2 l_f(√(n)x̅^e+θ̅) does not depend on N, W.

Similarly, we can provide a bound on g_t(u_t(0))-g_t(u_t^e).

    ∑_t=0^N-1 g_t(u_t(0))-g_t(u_t^e)   ≤∑_t=0^N-1∇ g_t(u_t(0)) u_t(0)-u_t^e
       ≤∑_t=0^N-1 l_g  u_t(0)-ξ_t u_t(0)-u_t^e
       ≤∑_t=0^N-1 l_g(√(n)u̅^e+ξ̅)  A(, :) x_t(0)-A(, :) x_t^e
       ≤∑_t=0^N-1 l_g(√(n)u̅^e+ξ̅)A(, :)x_t(0)-x_t^e 
       ≤ p^2l_g(√(n)u̅^e+ξ̅)A(, :)∑_t=0^N-1x_t-1^e-x_t^e

where the first inequality is by the convexity, the second inequality is because ξ_t is the global minimizer of g_t and g_t is l_g-smooth, the third inequality is by Assumption <ref>, Lemma <ref> and  Lemma <ref>, the fifth inequality is by  (<ref>). Notice that the constant term p^2l_g(√(n)u̅^e+ξ̅)A(, :) does not depend on N, W.

By (<ref>) and (<ref>), we complete the proof of the first inequality in the statement of Lemma <ref>:

    J(FOSS)-∑_t=0^N-1λ_t^e    ≤ c_1 ∑_t=0^N-1x_t-1^e-x_t^e + f_N(x_N(0))

where c_1 does not depend on N, W.



By defining x_N^e = θ_N, we can bound f_N(x_N(0)) by x_N(0)-x_N^e up to some constants because f_N(x_N(0))=f_N(x_N(0))-f_N(θ_N)  ≤l_f/2(√(n)x̅^e+ θ̅)  x_N(0)-x_N^e.  By the same argument as in (<ref>), we have x_N(0)-x_N^e = O(∑_t=0^N x_t-1^e- x_t^e ), where the big O hides some constant that does not depend on N, W. Consequently, 
 
    J(FOSS)-∑_t=0^N-1λ_t^e =O(  ∑_t=0^Nx_t-1^e-x_t^e)

 
 

 §.§ Proof of Lemma <ref>.

 The proof heavily relies on dynamic programming and the Bellman equations. For simplicity, we introduce a Bellman  operator  (f+g, h) defined by (f+g, h)(x)= min_u (f(x)+g(u)+h(Ax+Bu)). Now the Bellman equations can be written as (f+g, h^e)(x)=h^e(x)+λ^e for any x.
 
 
 We define a sequence of auxiliary functions S_k:  S_k(x)= h_k^e(x)+∑_t=k^N-1λ_t^e for k=0, …, N, where  h_N^e(x)=f_N(x). 
 
 We first provide a recursive equation for S_k. By Bellman equations, we have h_k^e(x)+λ_k^e=(f_k+g_k, h_k^e)(x) for 0≤ k ≤ N-1. Let π^e_k be the corresponding optimal control policy that solves the Bellman equations.  We have the following recursive relation for S_k  when 0≤ k ≤ N-1:
 
    S_k   (x) = (f_k+g_k, S_k+1-h^e_k+1+h_k^e)(x)

 where S_N(x)=f_N(x).
 
 Further, 	let V_k(x) denote the optimal cost-to-go function from k to N, then we  obtain a recursive equation for V_k by dynamic programming:
 
    V_k   (x) = (f_k+g_k, V_k+1)(x)= f_k(x)+g_k(π^*_k(x)) + V_k+1(Ax + B π^*_k(x))

	where 0≤ k ≤ N-1, and π^*_k denotes the optimal control policy and V_N(x)=f_N(x).
	
	Now, we are ready for a recursive inequality for S_k(x_k^*)-V_k(x_k^*). Let {x_k^*} denote the optimal trajectory, then x_k+1^* = Ax_k^*+Bπ^*_k(x_k^*). For any k=0,…, N-1,  
	
    S_k(x_k^*)-V_k(x_k^*)   = (f_k+g_k, S_k+1-h^e_k+1+h^e_k)(x_k^*)-(f_k+g_k, V_k+1)(x_k^*)
       ≤  f_k(x_k^*)+g_k(π^*_k(x_k^*)) + S_k+1(x_k+1^*)-h^e_k+1(x_k+1^*)+h_k^e(x_k+1^*)
         - (f_k(x_k^*)+g_k(π^*_k(x_k^*)) + V_k+1(x_k+1^*))
       =  S_k+1(x_k+1^*)-h^e_k+1(x_k+1^*)+h_k^e(x_k+1^*)-V_k+1(x_k+1^*)

	where the first inequality is because π_k^* is not optimal for the Bellman operator (f_k+g_k, S_k+1-h^e_k+1+h^e_k)(x_k^*).
	
	Summing over k=0, …, N-1 the recursive inequality for S_k(x_k^*)-V_k(x_k^*) yields
	
    S_0(x_0)-V_0(x_0)≤∑_k=0^N-1(h^e_k(x_k+1^*)-h^e_k+1(x_k+1^*))

	By subtracting h_0^e(x_0) on both sides,
	
    ∑_t=0^N-1λ^e_t-J^*≤∑_k=0^N-1(h^e_k(x_k+1^*)-h^e_k+1(x_k+1^*))-h^e_0(x_0)

	For the simplicity of notation, we define h^e_-1(x_0)=0 and x_0^*=x_0, then the bound can be written as
	
    ∑_t=0^N-1λ^e_t-J^*≤∑_k=0^N(h^e_k-1(x_k^*)-h^e_k(x_k^*))

 




 §.§ Proof of Lemma <ref>

The proof relies on the (strong) convexity and smoothness of the cost functions and the uniform upper bounds on θ_t, ξ_t.

	First of all, suppose there exists x̅^e such that x_t^e_2 ≤x̅^e for all 0≤ t ≤ N-1. We will bound u_t^e, x_t(0),u_t(0) by using x̅^e. Notice that the optimal steady state and the corresponding steady control satisfy: u_t^e= x_t^e,- A(, :)  x_t^e. If we can bound x_t^e by x_t^e ≤x̅^e for all t, u_t^e can be bounded accordingly:
	
    u_t^e   ≤ x_t^e, +  A(, :) x_t^e≤ x_t^e +  A(, :)x_t^e ≤ (1+ A(, :))x̅^e u̅^e
 
	Moreover, x_t(0) can also be bounded by x̅^e multiplied by some factors, because by Lemma <ref>, x_t(0)'s each entry is determined by some entry of x_s^e for s< t.  As  a result, for 0≤ t ≤ N
	
    x_t(0)_2 ≤√(n) x_t(0)_∞≤√(n)max_s< tx_s^e _∞≤√(n)max_s< tx_s^e _2≤√(n)x̅^e

	We can bound u_t(0) by  noticing that u_t(0)= x_t+1^(0)- A(, :) x_t(0) and 
	
    u_t(0)   ≤ x_t+1^(0) +  A(, :) x_t(0)≤ x_t+1(0) +  A(, :)x_t(0)
       ≤ (1+ A(, :))√(n)x̅^e= √(n)u̅^e


Next, it suffices to  prove x_t^e≤x̅^e for all t for some x̅^e. To prove this bound, we construct another (suboptimal) steady state: x̂_t= (θ^1_t, …, θ^1_t). Let û_t= x̂_t^-A(, :) x̂_t. It can be easily verified that (x̂_t, û_t) is indeed a steady state of the canonical-form system. Moreover, x̂_t and û_t can be bounded similarly as follows.

    x̂_t    ≤√(n) | θ_t^1| ≤√(n)θ_t _∞≤√(n)θ_t ≤√(n)θ̅
    û_t_2   ≤ (1+ A(, :)) x̂_t ≤ (1+ A(, :))√(n)θ̅


Now, we can bound x_t^e-θ_t.

    μ/2 x_t^e - θ_t ^2    ≤ f_t(x_t^e)- f_t(θ_t) + g_t(u_t^e)-g_t(ξ_t)  
       ≤ f_t(x̂_t)- f_t(θ_t) + g_t(û_t)-g_t(ξ_t)  
       ≤l_f/2x̂_t -θ_t^2 + l_g/2û_t- ξ_t^2 
       ≤ l_f (x̂_t ^2 + θ_t ^2) + l_g (û_t ^2 + ξ_t^2) 
       ≤ l_f (n θ̅^2 + θ̅^2) + l_g(( (1+ A(, :))√(n)θ̅)^2 + ξ̅)   c_5

where the first inequality is by f_t's strong convexity and g_t's convexity, the second inequality is because (x_t^e, u_t^e) is an optimal steady state, the third inequality is by the smoothness and ∇ f_t(θ_t)=∇ g_t(ξ_t)=0, the last inequality is by the bounds of x̂_t, û_t, θ_t, and ξ_t.




As a result, we have x_t^e - θ_t≤√( 2c_5/μ). Then,  we can bound x_t^e by x_t^e≤θ_t  + √( 2c_5/μ)≤θ̅+ √( 2c_5/μ)x̅^e for all t. It can be verified that x̅^e does not depend on N, W.
	
	



























§ LINEAR QUADRATIC TRACKING


In this section, we will provide a regret bound in Corollary <ref> for the general LQT  defined in Example <ref>. Based on this, we prove Corollary <ref>, which is a special case when Q_t, R_t are not changing.



 §.§ Regret bound on the general online LQT problems

Before the regret bound, we provide an important lemma to characterize the solution to the Bellman equations of the LQT problem.

One solution to the Bellman equations with stage cost 1/2(x-θ)^⊤ Q(x-θ)+1/2u^⊤ Ru can be represented by

    h^e(x)= 1/2(x-β^e)^⊤ P^e(x-β^e)

where P^e denotes the solution to the discrete-time algebraic Riccati equation (DARE) with respect to Q, R, A, B

    P^e=Q+A^⊤ (P^e-P^eB(B^⊤  P^eB+R)^-1B^⊤  P^e)A

and β^e= Fθ where F is a matrix determined by A, B, Q, R. 

The proof is in Appendix <ref>.

For simplicity of notation,  let P^e(Q, R) denote the solution to the DARE under the parameters Q, R, A, B and F(Q, R) denote the matrix in β^e=Fθ given parameters Q, R, A, B. Here we omit A, B in the arguments of the functions because they will not change in this paper.

In addition, we introduce the following useful notations: Q=μ_f I_n, Q̅= l_f I_n, R= μ_g I_m, R̅=l_g I_m for μ_f, μ_g >0, 0<l_f, l_g<+∞; and P̅=P^e(Q̅, R̅) and P=P^e(Q, R). Based on the notations above, we define some sets of matrices to be used later:
	
    ={Q|Q≤ Q ≤Q̅},
       = {R |R≤ R ≤R̅},
       ={P |P≤ P ≤P̅}.



Now, we are ready for the regret bound for the general LQT problem.

	Consider the LQT problem in Example <ref>. Suppose for t=0, 1, …, N-1, the cost matrices satisfy Q_t ∈, R_t ∈. Suppose the terminal cost function satisfies Q_N∈.[This additional condition is for technical simplicity and can be removed.]
	Then, the regret of RHTM with initialization FOSS can be bounded by
	
    (RHTM) = O( ζ^2(√(ζ)-1/√(ζ))^2K(∑_t=1^N( P^e_t-P^e_t-1 +  β_t^e-β^e_t-1) + ∑_t=0^N x_t-1^e- x_t^e ))

	where  K= (W-1)/p, x_-1^e=x_0, x_N^e=θ_N,   ζ is the condition number of the corresponding C( z), (x_t^e, u_t^e) is the  optimal steady state under cost Q_t, R_t, θ_t, P^e_t=P^e(Q_t, R_t) and β_t^e= F(Q_t, R_t)θ_t for t=0, …, N-1 and β_N^e=θ_N, P_N^e = Q_N.




Before the proof, we  introduce some supportive lemmas on the uniform bounds of P^e_t, β_t^e, x_t^* respectively. The intuition behind these uniform bounds is that the cost function coefficients Q_t, R_t, θ_t are all uniformly bounded by Assumption <ref> and <ref>.
The proofs are technical and deferred to Appendix <ref>.
	
		
	For  any Q_t ∈, R_t ∈, Q_N ∈, there exists x̅ that does not depend on t, N, W, such that 
	
    x_t^*_2 ≤x̅,     ∀ 0 ≤ t ≤ N.





  For any  Q∈, R∈, any θ≤θ̅, there exists β̅≥ 0 that does not depend on N and only depends on A, B, l_f, μ_f, l_g, μ_g, θ̅, such that max(θ̅,β^e )≤β̅, where β^e is defined in Lemma <ref>.
  
  
  For any  Q∈, R∈, we have P^e=P^e(Q, R) ∈. Consequently, P^e≤υ_max(P̅), where υ_max(P̅) denotes the largest eigenvalue of P̅.



Now, we are ready for the proof of Corollary <ref>. 

By Theorem <ref>, we only need to bound ∑_t=0^N(h^e_t-1(x_t^*)-h^e_t(x_t^*)). 
By definition, P_N^e=Q_N, β_N^e=θ_N, h_N^e(x)=f_N(x), so we can write h_t^e(x)= 1/2(x-β_t^e)^⊤  P_t^e(x-β_t^e) for 0 ≤ t ≤ N. 


For 0≤ t ≤ N-1, we split h_t^e(x_t+1^*)-h_t+1^e(x_t+1^*) into two parts.

    h_t^e(x_t+1^*)-h_t+1^e(x_t+1^*)   =1/2(x_t+1^*-β_t^e)^⊤ P_t^e(x_t+1^*-β_t^e)-1/2(x_t+1^*-β_t+1^e)^⊤ P_t+1^e(x_t+1^*-β_t+1^e)
       = 1/2(x_t+1^*-β_t^e)^⊤ P_t^e(x_t+1^*-β_t^e)-1/2(x_t+1^*-β_t+1^e)^⊤ P_t^e(x_t+1^*-β_t+1^e)_Part 1
        + 1/2(x_t+1^*-β_t+1^e)^⊤ P_t^e(x_t+1^*-β_t+1^e)-1/2(x_t+1^*-β_t+1^e)^⊤ P_t+1^e(x_t+1^*-β_t+1^e)_Part 2

Part 1 can be bounded by the following

    Part 1    = 1/2(x_t+1^*-β_t^e+  x_t+1^*-β_t+1^e )^⊤ P_t^e(x_t+1^*-β_t^e- (x_t+1^*-β_t+1^e))
       ≤1/2 x_t+1^*-β_t^e+  x_t+1^*-β_t+1^e _2P_t^e_2 β_t+1^e-β_t^e_2 
       ≤ (x̅ + β̅)υ_max(P̅) β_t+1^e-β_t^e_2

where the last inequality is by Lemma <ref>, <ref> <ref>.


Part 2 can be bounded by the following when 0≤ t ≤ N-1,

    Part 2    = 1/2(x_t+1^*-β_t+1^e)^⊤ (P_t^e-P_t+1^e)(x_t+1^*-β_t+1^e)
       ≤1/2x_t+1^*-β_t+1^e_2^2 P_t^e-P_t+1^e_2 ≤1/2(x̅+β̅)^2  P_t^e-P_t+1^e_2



Therefore, we have

    ∑_t=0^N(h^e_t-1(x_t^*)-h^e_t(x_t^*))   ≤∑_t=0^N-1(h_t^e(x_t+1^*)-h_t+1^e(x_t+1^*)) 
       =O( ∑_t=0^N-1( β_t+1^e-β_t^e_2+   P_t^e-P_t+1^e_2))

where the first inequality is by  h_0^e(x) ≥ 0 and h_-1^e(x)=0. 
Thus, by  Theorem <ref>, we have 

    Regret(RHTM) = O( ζ^2(√(ζ)-1/√(ζ))^2K(∑_t=1^N( P^e_t-P^e_t-1 +  β_t^e-β^e_t-1) + ∑_t=0^N x_t-1^e- x_t^e ))








 §.§ Proof of Corollary <ref>


Roughly speaking, the proof is mostly by applying Corollary <ref> and by showing β_t^e-β_t-1^e and x_t^e-x_t-1^e can be bounded by θ_t-θ_t-1 up to some constants and P_t^e-P^e_t-1=0 in the LQT problem (<ref>) where Q and R are not changing. However, directly applying  the results in Theorem <ref> and Corollary <ref> will result in some extra constant terms because some inequalities used to derive  the bounds in Theorem <ref> and Corollary <ref> are not necessary when Q, R are not changing. Therefore, we will need some intermediate results in the proofs of Theorem <ref> and Corollary <ref>  to prove Corollary <ref>.

Firstly, by  Lemma <ref> and Lemma <ref>, we have

    J(FOSS)- J^*     = J(FOSS)-∑_t=0^N-1λ_t^e + ∑_t=0^N-1λ_t^e-J^*
       ≤c_1 ∑_t=0^N-1 x_t-1^e- x_t^e _Part I+ ∑_t=0^N-1(h^e_t(x_t+1^*)-h^e_t+1(x_t+1^*))_Part II +f_N(x_N(0)) -h^e_0(x_0)_Part III

We are going to bound each part by ∑_t θ_t -θ_t-1 in the following. 


Part I: We will bound Part I  by ∑_t θ_t-θ_t-1 through showing that x_t^e=F_1F_2θ_t for some matrices F_1, F_2. The representation of x_t^e relies on Lemma <ref>.


By Lemma <ref>, any steady state (x,u) can be represented as a matrix multiplied by z: 

    x   =(z^1, …, z^1_p_1, z^2, …, z^2_p_2, …, z^m, …, z^m_p_m)^⊤ F_1z 
    
    u    = (z^1, …, z^m)^⊤ -A(, :) x  = (I_m-A(, :) F_1) z
 
where F_1∈^n,m is a binary matrix with full column rank. 

Consider cost function 1/2(x-θ)^⊤ Q(x-θ)+ 1/2u^⊤  R u. By the steady-state representation above, the optimal steady state can be solved by the following unconstrained optimization:

    min_z (F_1z-θ)^⊤ Q (F_1z-θ)+ z^⊤  (I-A(, :) F_1)^⊤ R(I-A(, :) F_1) z

Since F_1 is full column rank, the function is strongly convex and has the unique solution 

    z^e=F_2θ
  
where F_2 = (F_1^⊤  Q F_1+ (I-A(, :) F_1)^⊤ R(I-A(, :) F_1))^-1F_1^⊤ Q. Accordingly, the optimal steady state can be represented as  

    x^e= F_1 F_2θ,      u^e= (I_m-A(, :) F_1) F_2 θ.
 Consequently, when 1 ≤ t ≤ N-1,x_t^e-x_t-1^e≤F_1 F_2θ_t -θ_t-1.
When t=0, x_0^e -x_-1^e≤ F_1 F_2θ_0 -θ_-1 holds since x_-1^e= x_0=θ_-1=0.
Combining the upper bounds above, we have

    Part I =O( ∑_t=0^N-1θ_t -θ_t-1)





Part II: 
By (<ref>) in the proof of Corollary <ref>, and by noticing that P_t^e=P^e(Q, R) does not change, we have

    ∑_t=0^N-1(h_t^e(x_t+1^*)-h_t+1^e(x_t+1^*))    = O(∑_t=0^N-1β_t+1^e-β_t^e)

By Lemma <ref>, β_t^e = F(Q, R) θ_t  for 0≤ t ≤ N-1. In addition, since β_N^e=θ_N=0 as defined in (<ref>) and Corollary <ref>, we can also write  β_N^e=F(Q, R) θ_N.
Thus, 

    Part II=O(∑_t=0^N-1β_t+1^e-β_t^e) =O(  ∑_t=1^Nθ_t -θ_t-1)






Part III: 
By our condition for the terminal cost function, we have f_N(x_N(0))= 1/2(x_N(0)-β_N^e)^⊤  P^e (x_N(0)-β_N^e). By Lemma <ref>, we have h_0^e(x_0)= 1/2(x_0-β_0^e)^⊤  P^e (x_0-β_0^e). So Part III can be bounded by 

    Part III    = 1/2(x_N(0)-β_N^e)^⊤  P^e (x_N(0)-β_N^e)-1/2(x_0-β_0^e)^⊤  P^e (x_0-β_0^e)
        = 1/2(x_N(0)-β_N^e+ x_0-β_0^e)^⊤  P^e (x_N(0)-β_N^e-(x_0-β_0^e)) 
       ≤1/2 x_N(0)-β_N^e+ x_0-β_0^e P^e x_N(0)-β_N^e-(x_0-β_0^e)
       ≤1/2( √(n)x̅^e+ β̅+ β̅) P^e ( x_N(0)- x_0 + β_N^e -β_0^e)

where the last inequality is by x_0=0, Lemma <ref>, Lemma <ref>.

Next we will bound x_N(0)- x_0 and  β_N^e -β_0^e respectively. Firstly,  by β_t^e=F(Q,R)θ_t in Lemma <ref>, we have

    β_N^e -β_0^e   ≤∑_t=0^N-1β_t+1^e-β_t^e≤ F(Q,R)∑_t=0^N-1θ_t+1-θ_t

Secondly, we will bound x_N(0)- x_0. 

    x_N(0)- x_0   ≤ x_N(0)- x_N-1^e +  x_N-1^e- x_0
       ≤ x_N(0)- x_N-1^e + ∑_t=0^N-1x_t^e - x_t-1^e
       ≤ x_N(0)- x_N-1^e + F_1 F_2∑_t=0^N-1θ_t - θ_t-1

where the second inequality is by x_0^e=x_0, the third inequality is by (<ref>).


Next, we will focus on x_N(0)- x_N-1^e. By Lemma <ref>, 

    x_N(0)   = (z_N-p_1^e,1, …, z_N-1^e,1, z_N-p_2^e,2, …, z_N-1^e,2, …, z_N-p_m^e,m, …, z_N-1^e,m)^⊤
    
        x_N-1^e   = (z_N-1^e,1, …, z_N-1^e,1, z_N-1^e,2, …, z_N-1^e,2, …, z_N-1^e,m, …, z_N-1^e,m)^⊤

As a result,

    x_N(0)-x_N-1^e^2   ≤ z_N-2^e -z_N-1^e^2 + … +  z_N-p^e - z_N-1^e^2
        =F_2^2(θ_N-2-θ_N-1^2 + … + θ_N-p-θ_N-1^2)

where the equality is by (<ref>).
Taking square root on both sides yields

    x_N(0)-x_N-1^e   ≤F_1√(θ_N-2-θ_N-1^2 + … + θ_N-p-θ_N-1^2)
       ≤F_2 (θ_N-2-θ_N-1 + … + θ_N-p-θ_N-1)
       ≤F_2 (p-1) ∑_t=N-p^N-2θ_t+1-θ_t

Combining the bounds above leads to

    Part III=O(∑_t=0^N-1θ_t+1 -θ_t)

The proof is completed by summing up the bounds for Part I, II, III.





















§ PROOF OF THEOREM <REF>






Proof intuition:
By the problem transformation in Section <ref>, the fundamental limit of the online control problem is equivalent to the fundamental limit of the online convex optimization problem with objective C( z). Therefore, we will focus on C( z). Since the lower bound is for the worst case scenario, we only need to construct some tracking trajectories {θ_t} for Theorem <ref> to hold.   However, it is generally difficult to construct the tracking trajectories, so we consider randomly generated θ_t and show that the regret in expectation can be lower bounded. Then, there must exist some realization of the randomly generated {θ_t} such that the regret lower bound holds. 


 








Formal proof:





Step 1: construct LQ tracking. For simplicity, we construct a single-input system with n=p and  A ∈^n,n and B∈^n× 1 as follows: [It is easy to generalize the construction to multi-input case by constructing m decoupled subsystems. ]

    A= ([ 0 1 ⋯ 0; ⋮ ⋱ ⋱;     0 1; 1 0 ⋯ 0 ]),    B =(
    [ 0; ⋮; 0; 1 ])

(A, B) is controllable because (B, AB, …, A^p-1B) is full rank. A's controllability index is p=n. 

Next, we construct Q and R. For any ζ>1 and p, define δ = 4/(ζ-1)p. Let Q= δ I_n  and R= 1 for 0 ≤ t ≤ N-1. Let P^e=P^e(Q, R) be the solution to the DARE. The next lemma shows that P^e is a diagonal matrix and its diagonal entries can be characterized.
 
    Let P^e denote the solution to the DARE determined by A, B, Q, R defined above. Then P^e satisfies the form
    
    P^e= ([ q_1   0   ⋯   0;   0 q_2   ⋯   0;           ⋱;   0   ⋯ q_n ]),

where q_i = q_1+(i-1)δ for 1 ≤ i ≤ n and δ < q_1 < δ +1.


The DARE exists a unique positive definite solution <cit.>.
Suppose the solution is diagonal and substitute it in the DARE as follows.
 
    P^e
       = Q+A^⊤ (P^e-P^eB(B^⊤ P^eB+R)^-1B^⊤ P^e)A 
    ([ q_1   0   ⋯   0;   0 q_2   ⋯   0;           ⋱;   0   ⋯ q_n ])
              = ([ q_n/(1+q_n)+δ             0             ⋯             0;             0         q_1+δ             ⋯             0;                                         ⋱;             0             ⋯       q_n-1+δ ])

    So we have q_i = q_i-1+δ for 1 ≤ i ≤ n-1, and q_n/(1+q_n)+δ=q_1= q_n-(n-1)δ. Thus, q_n = nδ + √(n^2δ^2 + 4n δ)/2 >nδ. It is straightforward that q_1= q_n-(n-1)δ >δ >0, and q_1< δ +1 by  q_n/(1+q_n)<1. So we have found the unique positive definite solution to the DARE.



Next, we will construct θ_t. Let θ_0=θ_N=β_N^e=0 for simplicity. For θ_t when 1≤ t ≤ N-1, we divide the N-1 stages into E epochs, each with length Δ = N-1/L_N/2θ̅, possibly except  the last epoch. This is possible because 1≤Δ≤ N-1 by the conditions in Theorem <ref>. Thus, E= N-1/Δ. 
Let  be the first stage of the each epoch: = {1, Δ+1, …, (E-1)Δ +1 }. Let θ_t for t∈ independently and identically follow the  distribution below. 
    (θ_t^i=a)    =
    1/2    if  a=σ
    
    1/2    if  a=-σ,   i.i.d. for all i ∈[n], t∈,

where σ= θ̅/√(n). It can be easily verified that θ= θ̅ for any realization of this distribution, so Assumption <ref> is satisfied. Let the other  θ_t in each epoch be equal to the θ at the start of their corresponding epochs, i.e. θ_kΔ +1=θ_kΔ +2=⋯ =θ_(k+1)Δ, when k ≤ E-1, and θ_kΔ +1=… =θ_N-1 when k=E. The following inequalities show that the constructed {θ_t} satisfies the variation budget:

    ∑_t=0^N θ_t -θ_t-1   =θ_1-θ_0 + ∑_k=1^E-1θ_kΔ+1-θ_kΔ + θ_N-1-θ_N
       ≤θ̅+ 2(E-1)θ̅+ θ̅=2θ̅E
       ≤ 2θ̅L_N/2θ̅≤ 2θ̅L_N/2θ̅= L_N

where the first equality is 
by θ_0=θ_-1=θ_N=0, the first inequality is by θ_t =θ̅ when 1 ≤ t ≤ N-1, the second inequality is by Δ = N-1/L_N/2θ̅≥N-1/L_N/2θ̅,   and thus L_N/2θ̅≥N-1/Δ=E.


















The  total cost of our constructed LQ tracking problem is

    J( x, u)= ∑_t=0^N-1(δ/2x_t-θ_t^2 + 1/2u_t^2)+ 1/2x_N^⊤  P^e x_N



We will verify that C( z)'s condition number is ζ in Step 2.




Step 2: problem transformation and  the optimal solution z^*. By the problem transformation in Section <ref>, we let z_t =x_t^n, and the equivalent cost function C( z) is given below.

    C( z)=∑_t=0^N-1(δ/2∑_i=1^n ( z_t-n+i- θ_t^i)^2 +1/2(z_t+1- z_t-n+1)^2) +  1/2∑_i=1^n q_iz_N-n+i^2

and z_t =0 and θ_t=0 for t ≤ 0. 

Since C( z) is strongly convex, min C( z) admits a unique optimal solution, denoted as z^*, which is determined by the first-order optimality condition: ∇ C( z^*)=0. In addition, our constructed C( z) is a quadratic function, so there exists a matrix H ∈^N× N and a vector η∈^N such that ∇ C( z^*)=H z^*-η=0. By the partial gradients of C( z) below,

    ∂ C/∂ z_t= δ ( z_t - θ_t^n + z_t -θ_t+1^n-1 +… + z_t - θ_t+n-1^1) + z_t - z_t+n+z_t - z_t-n,  1 ≤ t ≤ N-n
       ∂ C/∂ z_t= δ ( z_t - θ_t^n +… + z_t - θ_N-1^n+t-N+1 ) + q_n+t-Nz_t+z_t-z_t-n,      N-n+1 ≤ t ≤ N

For simplicity and without loss of generality, we assume  that N/n is an integer. Then, by Lemma <ref>,
 H can be represented as  the block matrix below 

    H = (
    [ (δ n+2)I_n       -I_n          ⋯;       -I_n (δ n+2)I_n          ⋱           ;                     ⋱          ⋱       -I_n;                             -I_n (q_n+1)I_n ])∈^N × N.

 η is a linear combination of θ: for 1 ≤ t ≤ N, we have 
η_t = δ(  θ_t^n+… +   θ_t+n-1^1) = δ(e_n^⊤θ_t  + … + e_1^⊤θ_t+n-1) 
where e_1, …, e_n∈^n are standard basis vectors and θ_t=0 for t≥ N.





By Gergoskin's Disc Theorem and Lemma <ref>, H's condition number is (δ n+4)/δ n =ζ by our choice of δ in Step 1 and p=n. Thus we have shown that C( z)'s condition number is ζ.

Since H is strictly diagonally dominant with positive diagonal entries and nonpositive off-diagonal entries, H is invertible and its inverse, denoted by Y, is nonnegative. Consequently,  the optimal solution can be represented as z^*=Yη. Since η is linear in {θ_t }, z_t^* is also linear in {θ_t } and can be characterized by the following.

    z_t+1^*    = ∑_i=1^N Y_t+1,iη_i =δ∑_i=1^N Y_t+1,i∑_j=0^n-1 e_n-j^⊤θ_i+j
        = δ∑_k=1^N-1( ∑_i=1^n Y_t+1,i+k-ne_i^⊤) θ_k
       δ∑_k=1^N-1 v_t+1,kθ_k

where  θ_t=0 for t≥ N, Y_t+1,i=0 for i≤ 0, and v_t+1,k∑_i=1^n Y_t+1,i+k-ne_i^⊤.
 
 

In addition, we are able to show in the next lemma that Y has decaying row entries starting at the diagonal entries. The proof is technical and deferred to the Appendix F.1. 

When N/n is an integer, 
the inverse of H, denoted by Y, can be represented as a block matrix

    Y = (
    [     y_1,1 I_n     y_1,2 I_n             ⋯   y_1,N/n I_n;     y_2,1 I_n     y_2,2 I_n             ⋯   y_2,N/n I_n;             ⋮             ⋱             ⋱             ⋮;   y_N/n,1 I_n    y_N/n,2I_n             ⋯ y_N/n,N/n I_n ])

where y_t, t+τ≥1-ρ/δ n +2ρ^τ>0 for τ≥ 0 and ρ= √(ζ)-1/√(ζ)+1.






Step 3: characterize z_t+1(^z).
For any online control algorithm , we can define an equivalent online algorithm for z, denoted as ^z. ^z, at each time t, outputs z_t+1(^z)  based on the predictions and the history, i.e.,

    z_t+1(^z) = ^z( {θ_s }_s=0^t+W-1),    t≥ 0

For simplicity, we consider online deterministic algorithm.[The proof can be easily generalized to random algorithms] Notice that z_t+1 is a random variable because θ_1, …, θ_t+W-1 are random. Based on this observation and Lemma <ref>, we are able to provide a regret lower bound in Step 4.





Step 4: prove the regret lower bound on .
Roughly speaking, the regret occurs when something unexpected happens beyond the prediction window, that is,  at each t, the prediction window goes as far as t+W-1, but if θ_t+W changes from θ_t+W-1, the online algorithm cannot prepare for it, resulting in poor control and positive regret. 
By our construction, when t+W∈, θ_t+W changes from θ_t+W-1. To study such t, we  define a set  _1 = { 0 ≤ t ≤ N-W-1| t+W ∈}.
It can be shown that the cardinality of _1 can be lower bounded by L_N up to some constants:

    |_1| ≥1/18 θ̅L_N

The proof of (<ref>) is provided below.

    |_1|    = | { W≤ t ≤ N-1| t∈}|
       =| |-| { 1≤ t ≤ W-1| t∈}|
        = N-1/Δ- W-1/Δ
       ≥N-W/Δ
       ≥1/2N-W/Δ
       ≥1/2N-W/N-1 + L_N/2θ̅L_N/2θ̅
       ≥1/2N-1/3N/N-1+N+1/2L_N/2θ̅≥1/6L_N/2θ̅
       ≥1/62/3L_N/2θ̅= 1/18L_N/θ̅

where the first inequality is by the definition of the ceiling and floor operators, the second inequality is by N-W/Δ≥ 1 under the conditions on N, W,  L_N in Theorem <ref>, the third inequality is by Δ = N-1/L_N/2θ̅≤N-1/L_N/2θ̅+1, the fourth inequality is by L_N ≤ (2N+1)θ̅ in Theorem 3's statement, the last inequality is by L_N ≥ 4 θ̅ in Theorem 3's statement.












Moreover, we can show  in Lemma <ref> that, for all  t ∈_1,  the online decision z_t+1(^z) is different from the optimal solution  z_t+1^* and the difference is lower bounded,   






 For any online algorithm ^z,
    when t ∈_1, 
    
    |z_t+1(^z)-z_t+1^*|^2 ≥ c_10σ^2ρ^2K

    where c_10 is a constant determined by A, B, n, Q, R constructed above and ρ= √(ζ)-1/√(ζ)+1.

The proof is provided in Appendix F.2.

The lower bound on the difference between the online decision and the optimal decision results in a lower bound on the regret. 
By the nδ-strong convexity of C( z),

    (C( z(^z))-C( z^*) )   ≥δ n/2∑_t∈_1 |z_t+1(^z)-z_t+1^*|^2
       ≥ |_1| c_10σ^2ρ^2K
       ≥L_N/18 θ̅ c_10σ^2ρ^2K = Ω(L_N ρ^2K)


By the equivalence between  and ^z, we have
J()- J^* = Ω(ρ^2KL_N). 
By the property of expectation, there must exist some realization of the random {θ_t} such that J()- J^* = Ω(ρ^2KL_N), where ρ= √(ζ)-1/√(ζ)+1. This completes the proof. 













 §.§ Proof of Lemma <ref>



	
	Since H is a block matrix
	
    H = (
    	[ (δ n+2)I_n       -I_n          ⋯;       -I_n (δ n+2)I_n          ⋱           ;                     ⋱          ⋱       -I_n;                             -I_n (q_n+1)I_n ])

	its inverse matrix Y can also be represented as a block matrix. Moreover, let 
	
    H_1 = ([ δ n +2     -1      ⋯      0;     -1 δ n +2      ⋱      0;      ⋮      ⋱      ⋱      ⋮;      0      ⋯     -1  q_n+1 ])
 and define Y̅=(H_1)^-1=(y_ij)_i,j =1^N/n. Then the inverse matrix Y can be represented as the block matrix: Y=(y_ijI_n)_i,j =1^N/n. 
	
	Now, it suffices to provide a lower bound on y_ij.



	Since H_1 is a symmetric positive definite tridiagonal matrix, by <cit.>, the inverse has an explicit formula given by (H_1)^-1_ij= a_ib_j and 
	
    a_t    = ρ/1- ρ^2( 1/ρ^t   -ρ^t)
    
    	b_t    = c_3 1/ρ^N-t + c_4 ρ^N-t
    
    	c_3     = b_N( (q_n+1)ρ -ρ^2/1-ρ^2)
    
    	c_4     = b_N 1-(q_n+1)ρ/1-ρ^2
    
    	b_N    = 1/- a_N-1+ (q_n+1)  a_N

	
	In the following, we will show y_t,t+τ= a_t b_t+τ≥1-ρ/δ n +2ρ^τ when τ≥ 0.  Firstly, it is easy to verify that
	
    ρ^t   a_t = ρ/1-ρ^2(1-ρ^2t)≥ρ

	since t≥ 1 and ρ<1.
	
	Secondly, we bound b_N in the following way:
	
    ρ^-N b_N  = 1/(q_n+1)(1-ρ^2N)-(ρ-ρ^2N-1) 1-ρ^2/ρ≥1/(δ n+2) 1-ρ^2/ρ

	because 0<(q_n+1)(1-ρ^2N)-(ρ-ρ^2N-1) ≤ (δ n+2) by nδ <q_n <nδ +1 in Lemma <ref>.
	
	
	Thirdly, we bound b_t+τ. When 1-(q_n+1)ρ≥ 0
	
    ρ^N-t-τb_t+τ   =  b_N( (q_n+1)ρ -ρ^2/1-ρ^2) +  b_N 1-(q_n+1)ρ/1-ρ^2ρ^2(N-t-τ)
       ≥ b_N( (q_n+1)ρ -ρ^2/1-ρ^2) 
       ≥ b_N( (δ n+1)ρ -ρ^2/1-ρ^2) 
        = 1-ρ/1-ρ^2b_N

	where the first inequality is by
	 1-(q_n+1)ρ≥ 0, the second inequality is by qn >nδ in Lemma <ref>, and
	the last equality is by ρ^2 - (δ n+2)ρ +1=0.
	
	When  1-(q_n+1)ρ < 0
	
    ρ^N-t-τb_t+τ   =  b_N( (q_n+1)ρ -ρ^2/1-ρ^2) +  b_N 1-(q_n+1)ρ/1-ρ^2ρ^2(N-t-τ)
       ≥ b_N( (q_n+1)ρ -ρ^2/1-ρ^2) +  b_N 1-(q_n+1)ρ/1-ρ^2
       ≥ b_N ≥1-ρ/1-ρ^2b_N

	where the first inequality is by 1-(q_n+1)ρ < 0, ρ≤ 1, the second inequality is by ρ^2(N-t-τ)≤ 1. Thus, we obtained a lower bound for b_t+τ.
	
	Combining bounds of a_t, b_t+τ, b_N together yields
	
    y_t,t+τ= a_tb_t+τ   ≥ρ b_N1-ρ/1-ρ^2ρ^τ-N≥1-ρ/(δ n+2)ρ^τ

	





 §.§ Proof of Lemma <ref>




    By our construction, θ_t is random,  z_t+1^ is also random and its randomness is provided by θ_1, …, θ_t+W-1, while z_t+1^* is determined by all θ_t. When t∈_1, 
	
    |z_t+1^-z^*_t+1|^2    =  | z_t+1^- δ∑_i=1^N-1 v_t+1,iθ_i|^2  
        =  |z_t+1^- δ∑_i=1^t+W-1v_t+1,iθ_i ^2 +δ^2  | ∑_i=t+W^N-1 v_t+1,iθ_i|^2 
       ≥δ^2  | ∑_i=t+W^N-1 v_t+1,iθ_i|^2,

	where the first equality is by (<ref>), the second equality is by θ_τ=0 for all τ, and θ_t+W, …, θ_N are independent of θ_1, …, θ_t+W-1 when t∈_1.
	
	Further,
	
    | ∑_i=t+W^N-1 v_t+1,iθ_i|^2    =  | ∑_i=t+W^t+W+Δ-1 v_t+1,iθ_t+W|^2+…+ | ∑_i=(E-1)Δ+1^N-1 v_t+1,iθ_(E-1)Δ+1|^2
        = ∑_i=t+W^t+W+Δ-1 v_t+1,i^2 σ^2+… + ∑_i=(E-1)Δ+1^N-1 v_t+1,i^2 σ^2
       ≥σ^2 ∑_i=t+W^N-1v_t+1,i^2  
       = σ^2 ∑_i=t+W^N-1 (∑_k=0^n-1 Y_t+1,i-k^2)  ≥σ^2 ∑_i=t+1+W-n^N-1 Y_t+1,i^2  
       =σ^2 ∑_i=t+1+W-n^N Y_t+1,i^2

	where the first equality is because the theta in one epoch are equal by our construction,
	the second equality is because cov(θ_τ)=σ^2 I_n, 
	the first inequality is because  the entries of v_t+1,i  are nonnegative, the third equality is by the definition of v_t+1,i in (<ref>), and
	the last equality is because when t ∈_1, Y_t+1, N=0.
	
		When 1≤ W≤ n, ∑_i=t+1+W-n^N Y_t+1,i^2   ≥ Y_t+1,t+1^2= Y_t+1, t+1+n W-1/n. When W>n, ∑_i=t+1+W-n^N Y_t+1,i^2  ≥ Y_t+1,t+1+nW-n/n^2. Moreover, when W≥ 1, W-n/n=W-1/n. In summary, for  W≥ 1, 
	
    ∑_i=t+1+W-n^N Y_t+1,i^2   ≥  Y_t+1,t+1+nW-1/n^2 ≥ρ^2 K(1-ρ/δ n +2)^2

	where the last inequality is by Lemma <ref>. This completes the proof. 








§ PROOFS OF THE LQT'S PROPERTIES USED IN APPENDIX <REF>


In this section, we provide proofs for the properties of LQ tracking (LQT) used in Appendix <ref>. 









 §.§ Preliminaries: dynamic programming for finite-horizon LQT


In this section, we consider a discrete time LQ tracking problem with time-varying cost functions and time-invariant dynamical system:

    min_x_t,u_t     1/2∑_t=0^N-1[ (x_t-θ_t)^⊤  Q_t (x_t-θ_t)+ u_t^⊤  R_t u_t ] + 1/2 (x_N-θ_N)^⊤  Q_N(x_N-θ_N)
    s.t.      x_t+1=Ax_t + Bu_t,      t=0, …, N-1

where x_0=0 for simplicity. 

The problem can be solved by dynamic programming.


	Consider a finite-horizon time-varying LQ tracking problem. Let V_t(x_t) be the cost to go from k=t to k=N, then 
	
    V_t(x_t)= 1/2 (x_t-β_t)^⊤  P_t(x_t-β_t)+1/2∑_k=t^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)

	for t=0, …, N. The parameters can be obtained by
	
    P_t    = Q_t + A^⊤  M_t A,    t=0, …, N-1,    P_N=Q_N
    
    	 M_t    = P_t+1-P_t+1B(R_t+B^⊤  P_t+1B)^-1B^TP_t+1,    t=0, …, N-1
    β_t     = (Q_t + A^⊤  M_t A)^-1 (Q_t θ_t + A^⊤  M_t β_t+1),    t=0, …, N-1
    β_N    = θ_N
    
    H_t	     = M_t - M_t A(Q_t +A^⊤ M_t A)^-1 A^⊤  M_t ,    t=0, …, N-1

	
	The optimal controller is 
	
    u_t^* = - K_t x_t + K_t' β_t+1,    t=0, …, N-1

	where the parameters are 
	
    K_t=(R_t + B^⊤  P_t+1B)^-1 B^⊤  P_t+1A
       K_t'=(R_t + B^⊤  P_t+1B)^-1 B^⊤  P_t+1

	
	There is another way to write the optimal controller: 
	
    u_t^* =- K_t x_t + K_t^αα_t+1   t=0, …, N-1

	where the parameters are 
	
    K_t^α   = (R_t + B^⊤  P_t+1B)^-1 B^⊤
    α_t    = P_t β_t
    α_t    = Q_t θ_t +  (A-BK_t)^⊤α_t+1,    t=0, …, N-1
    α_N   =P_Nθ_N

	
	



    The proof is straightforward by following dynamic programming procedures.
    
    Firstly, it is direct to verify that V_N(x_N)= 1/2(x_N-θ_N)^⊤ Q_N(x_N-θ_N).
    Then, suppose the claim of Theorem <ref> is true at t+1, we will verify the stage t in the following.
    
    V_t(x_t)=   min_u_t[1/2(x_t-θ_t)^⊤ Q_t (x_t-θ_t) +1/2 u_t^⊤ R_t u_t+ V_t+1(Ax_t+Bu_t)]
    
             =   1/2min_u_t[(x_t-θ_t)^⊤ Q_t (x_t-θ_t) + u_t^⊤ R_t u_t + (Ax_t+Bu_t-β_t+1)^⊤ P_t+1(Ax_t+Bu_t -β_t+1)
         +  ∑_k=t+1^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)
            ]
    
             =    1/2(Ax_t-β_t+1)^⊤ (P_t+1-P_t+1B(R_t +B^⊤ P_t+1B)^-1B^⊤ P_t+1)(Ax_t-β_t+1)
         + 1/2(x_t-θ_t)^⊤ Q_t (x_t-θ_t)+ 1/2∑_k=t+1^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)
    
            =   1/2(Ax_t-β_t+1)^⊤ M_t (Ax_t-β_t+1)  + 1/2(x_t-θ_t)^⊤ Q_t (x_t-θ_t)
       + 1/2∑_k=t+1^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)
    
            =   1/2(x_t-β_t)^⊤ P_t(x_t-β_t) - 1/2(Q_t θ_t+ A^⊤ M_t β_t+1)^⊤(Q_t +A^⊤ M_t A)^-1(Q_t θ_t+ A^⊤ M_t β_t+1)
       + 1/2θ_t^⊤ Q_t θ_t + 1/2β_t+1^⊤ M_t β_t+1+
            1/2∑_k=t+1^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)
    
            =   1/2(x_t-β_t)^⊤ P_t(x_t-β_t)+
            1/2∑_k=t^N-1  (Aθ_k  - β_k+1)^⊤  H_k (Aθ_k  - β_k+1)

    where the third equality is by noticing that the optimal control input is
    
    u_t^*= -(R_t +B^⊤ P_t+1B)^-1B^⊤ P_t+1(Ax_t-β_t+1)=-K_tx_t+K_t'β_t+1,

    the fourth equality is by M_t's definition, the fifth equality is by combining the two quadratic terms of x_t as one quadratic term with a constant, and the last equality is by definition.


























 §.§ Proof of Lemma <ref>

 In the following, we first prove that the recursive solution P_t to the finite-horizon LQT is bounded. Then,  we can prove Lemma <ref> by taking limits.

Consider a finite-horizon time-varying LQT problem.	For any N, any 0≤ t ≤ N, any Q_t ∈, R_t ∈, Q_N ∈, we have P_t ∈ where P_t is defined in Theorem <ref>.


In the following, we use the notations and definitions introduced in Appendix <ref>.1 and Theorem <ref>.

Since P_t does not depend on θ_t, we let θ_t=0 and consider the LQR problem for simplicity.
	Since Q≤ Q_t ≤Q̅, R≤ R_t ≤R̅, for 0≤ t ≤ N-1 and P≤ Q_N ≤P̅, we have for any x_t, u_t, k, Q_t, R_t, Q_N,
	
    ∑_t=k^N-1 (x_t^⊤  Q_t x_t +u_t^⊤  R_t u_t)+ x_N^⊤ Q_Nx_N ≤∑_t=k^N-1 (x_t^⊤Q̅ x_t +u_t^⊤R̅ u_t)+ x_N^⊤P̅x_N
       ∑_t=k^N-1 (x_t^⊤  Q_t x_t +u_t^⊤  R_t u_t)+ x_N^⊤ Q_Nx_N ≥∑_t=k^N-1 (x_t^⊤Q x_t +u_t^⊤R u_t)+ x_N^⊤Px_N

	Taking minimum over all feasible trajectories on both sides yields
	
    min∑_t=k^N-1 (x_t^⊤  Q_t x_t +u_t^⊤  R_t u_t)+ x_N^⊤ Q_Nx_N ≤min∑_t=k^N-1 (x_t^⊤Q̅ x_t +u_t^⊤R̅ u_t)+ x_N^⊤P̅x_N
       min∑_t=k^N-1 (x_t^⊤  Q_t x_t +u_t^⊤  R_t u_t)+ x_N^⊤ Q_Nx_N ≥min∑_t=k^N-1 (x_t^⊤Q x_t +u_t^⊤R u_t)+ x_N^⊤Px_N

	Notice that the left-hand-side terms of both inequalities are equal to x_k^⊤ P_k x_k. Moreover, notice that
	
    x_k^⊤P̅ x_k=min_x_t+1=Ax_t+Bu_t∑_t=k^N-1 (x_t^⊤Q̅ x_t +u_t^⊤R̅ u_t)+ x_N^⊤P̅x_N

	because P̅=P^e(Q̅, R̅) is the solution to the DARE. The same holds for P.
	Therefore, we have
	
    x_k^⊤P x_k ≤ x_k^⊤  P_k x_k≤ x_k^⊤P̅ x_k

	for any x_k. Thus, P≤ P_k ≤P̅, i.e. P_k ∈.



 In the following, we use the notations and definitions introduced in Appendix <ref>.1 and Theorem <ref>. 
Since P^e is not influenced by θ_t, we let θ_t=0 for simplicity. Consider a finite-horizon LQR problem: ∑_k=0^N-1 (x_k^⊤  Q x_k +u_k^⊤  R u_k)+ x_N^⊤ Q_N x_N, where Q_N∈. By Lemma <ref>, we have P_k ∈. Since P_k → P^e as k→ -∞ <cit.>, and since  is a closed set <cit.>,  we have P^e ∈. Since P^e and P̅ are positive definite, we have P^e _2≤υ_max(P̅). 



	




























 §.§ Proof of Lemma <ref>

In the following, we will provide and prove an enhanced version of Lemma <ref> with detailed characterization of the solution to the Bellman equations in Proposition <ref>. 



	Suppose (A, B) is controllable, Q,R are positive definite. The optimal average cost λ^e does not depend on the initial state x_0 and is equal to
	
    λ^e =  1/2(Aθ-β^e)^⊤  H^e(Aθ-β^e),

	where M^e=P^e-P^eB(R+B^⊤ P^eB)^-1B^⊤ P^e and H^e=M^e-M^eA(Q+A^⊤ M^eA)^-1A^⊤ M^e.
	
	In addition, a bias function of the Bellman equations h^e(x)+λ^e= min_u (f(x)+g(u)+h^e(Ax+Bu)) can be represented by
	
    h^e(x)=1/2(x-β^e)^⊤  P^e(x-β^e).

	where P^e=P^e(Q,R).
	
	The optimal controller is 
    u=-K^e x +K'β^e

	 where 
	 K^e=(R+B^⊤ P^eB)^-1B^⊤ P^eA, K'=(R+B^⊤ P^eB)^-1B^⊤ P^e, and β^e satisfies
	 
    β^e=(P^e)^-1α^e=Fθ

	 where α^e= Qθ+(A-BK^e)^⊤α^e and thus 
	 F= (P^e)^-1(I-(A-BK^e)^⊤ )^-1Q. 
	 
	 




It is easy to see that the formulas of λ^e, h^e(x), and the optimal controller are  the limits of the corresponding formulas or the limiting solutions to the corresponding iterative equations  under fixed Q, R, θ in Theorem <ref>. However, to formally prove these formulas, we still need to prove the existence of the limits, which is the focus of the following proof. In particular, the proof consists of three parts: i) verify the formula of the optimal average cost λ^e, ii) verify the formula of the  bias function h^e(x), iii) verify the formula of the optimal controller.

Part i): Verify the formula of λ^e. 
		Consider a finite horizon LQT problem:
	
    min_x_t,u_t     1/2∑_t=0^N-1[ (x_t-θ)^⊤  Q(x_t-θ)+ u_t^⊤  R u_t ] 
    s.t.      x_t+1=Ax_t + Bu_t,      t=0, …, N-1

	Given an initial state x_0, by Theorem <ref>, the total optimal cost in N time steps is 
	
    J_N^*(x_0) =1/2 (x_0-β_0)^⊤  P_0(x_0-β_0)+1/2∑_k=0^N-1  (Aθ  - β_k+1)^⊤  H_k (Aθ  - β_k+1)

	If we can show that  β_k →β^e and P_k→ P^e and H_k→ H^e as k → -∞, then, consequently, we will have 1/2 (Aθ  - β_k+1)^⊤  H_k (Aθ  - β_k+1) →1/2  (Aθ  - β^e)^⊤  H^e (Aθ  - β^e) as k→ -∞, and  bounded 1/2 (x_0-β_0)^⊤  P_0(x_0-β_0) for fixed x_0.  Then the formula of the optimal average cost in infinite horizon can be proved by   
	
    λ^e    =  lim_N → +∞[1/Nmin_x_t+1=Ax_t + Bu_t( 1/2∑_t=0^N-1( (x_t-θ)^⊤  Q(x_t-θ)+ u_t^⊤  R u_t ) ) ]
       =
    	lim_N → +∞[1/N(1/2 (x_0-β_0)^⊤  P_0(x_0-β_0)+1/2∑_k=0^N-1  (Aθ  - β_k+1)^⊤  H_k (Aθ  - β_k+1))]
       =1/2(Aθ-β^e)^⊤  H^e(Aθ-β^e),

		Therefore, it suffices to prove β_k →β^e, P_k→ P^e and H_k→ H^e as k → -∞. 
		
		By Proposition 4.4.1 <cit.>, P_k → P^e as k → -∞. Then, M_k → M^e  as k → -∞ since M_k is a continuous function of P_k by noticing that the matrix inverse operator is continuous when the matrix is invertible. Similarly, H_k → H^e as k → -∞ since H_k is a continuous function of M_k. In addition, K_k→ K^e, and K_k^α→ K^α and K'_k → K' as k → -∞ since K_k, K_k^α, K'_k are continuous functions of P_k. 
		
		To show β_k →β^e, we only need to show 
		α_k →α^e as k → -∞ since β_k = P_k^-1α_k. α_k satisfies the recursive equation α_k= Qθ +(A-BK_k+1)^⊤α_k+1. Since (A-BK_k)^⊤→ (A-BK^e)^⊤ as k → -∞ and (A-BK^e)^⊤ is a stable matrix, by the claim below, we can show α_k →α^e as k → -∞. Then, the proof of Part i) is completed.
		
		Claim: 
		If A_t→ A and A is stable, then the state of the system x_t+1=A_t x_t +η  will converge to x^s, where x^s=Ax^s+η, for any bounded initial value x_0.
		
The proof of the claim lemma is provided at the end of this subsection.


Part ii): Verify h^e(x)'s formula.
The proof is by showing the Bellman equations hold under the formulas of h^e(x) and λ^e provided in the statement of the lemma. 

    min_u    [1/2 (x-θ)^⊤  Q(x-θ)+ 1/2u^⊤  R u + 1/2(Ax+Bu-β^e)^⊤  P^e(Ax+Bu-β^e)]
    
      =     1/2 (x-θ)^⊤  Q(x-θ)+ 1/2 (Ax-β^e)^⊤ M^e (Ax-β^e)
       +min_u 1/2(u+ K^e x-K'β^e)^⊤  (R+B^⊤ P^e B)( u +K^e x-K'β^e)
    
     =     1/2 (x-θ)^⊤  Q(x-θ) + 1/2(Ax -β^e)^⊤  M^e(Ax -β^e) 
    
     =     1/2(Aθ-β^e)^⊤  H^e(Aθ-β)+ 1/2(x-β^e)^⊤  P^e(x-β^e)

where the last equality is by Q+A^⊤ M^e A= P^e,  β^e = (P^e)^-1α^e,  α^e = Qθ +(A-BK^e)^⊤α^e and the formulas of K^e, M^e.


	


Part iii): Verify the formula of the optimal controller. We prove u=-K^e x +K'β^e is the optimal controller  by showing that the average cost by implementing this controller is no more than the optimal average cost λ^e. Let x_t,u_t be the state and control at t by implementing the controller u=-K^e x +K'β^e.
	
    1/N( 1/2∑_t=0^N-1[ (x_t-θ)^⊤  Q(x_t-θ)+ u_t^⊤  R u_t ]) 
       ≤1/N(1/2∑_t=0^N-1[ (x_t-θ)^⊤  Q(x_t-θ)+ u_t^⊤  R u_t ]+ 1/2(x_N-β^e)^⊤ P^e(x_N-β^e))
       = 1/N(1/2 (x_0-β^e)^⊤  P^e(x_0-β^e)+1/2∑_k=0^N-1  (Aθ  - β^e)^⊤  H^e (Aθ  - β^e))

where the last equality is by Theorem <ref>.
Taking N→ +∞ on both sides, we have

    lim_N→ +∞1/N1/2∑_t=0^N-1[ (x_t-θ)^⊤  Q(x_t-θ)+ u_t^⊤  R u_t ] 
    ≤1/2  (Aθ  - β^e)^⊤  H^e (Aθ  - β^e)=λ^e

This completes the proof.


Proof of Claim:
Define the error term d_t =x_t-x^s. The dynamics of d_t is d_t+1=Ad_t+w_t, where w_t=(A_t-A)(d_t+x^s). It suffices to show d_t → 0 as t→ +∞. In the following, we will first prove two facts, based on which we prove d_t → 0.

Fact 1: Consider a stable matrix A and a sequence of uniformly bounded vectors: v_t _2 ≤ D for all t. There exists a constant c_6>0 determined by A, such that, for any k=1, 2, …, 

    ∑_t=0^k-1 A^t v_t _2 ≤ c_6 D.


Proof of Fact 1: This is a consequence of the fact that exponential stability implies bounded-input-bounded-output stability. To see this, 	consider a system x_t+1=Ax_t+u_t with x_0=0. Since A is stable, the system is exponentially stable. By Theorem 9.4 <cit.>, the exponential stability implies the bounded-input-bounded-output stability. Thus, there exists c_6 such that x_k_2 ≤ c_6 D for any k and  any  input sequence satisfying u_t _2 ≤ D for all t. 

For any k≥ 0, consider inputs u_t=v_k-1-t for 0 ≤ t ≤ k-1, then x_k= ∑_t=0^k-1 A^tv_t. Since u_t ≤ D,  we have
x_k≤ c_6 D, which completes the proof.





Fact 2: There exists a constant D>0, such that max_t≥ 0(x^s, d_t )≤ D.

Proof of Fact 2:
Since A_t→ A, for  ϵ_1= 1/(4c_6), there exists N_1, such that when t≥ N_1, A_t-A≤ϵ_1. 
Since A is stable, we have A^t  → 0, so for  ϵ_2=1/2, there exists N_2, such that when t>N_2, A^t ≤ϵ_2. 
	
	
	 Let D=max(d_0, …, d_N_1+N_2, x^s). By definition, d_t≤ D for t≤ N_1+N_2. We can show  d_N_1+N_2+1≤ D in the following.
	
    d_N_1+N_2+1   = A^N_2+1d_N_1+ w_N_1+N_2 +Aw_N_1+N_2-1+… + A^N_2w_N_1
       ≤A^N_2+1_2 D +  w_N_1+N_2 +Aw_N_1+N_2-1+… + A^N_2w_N_1
       ≤ϵ_2 D + c_6 max_N_1≤ k ≤ N_1+N_2w_k
       ≤ϵ_2 D + 2c_6 ϵ_1 D  = (1/2+1/2)D=D

	where the second inequality is by Fact 1 and the definitions of ϵ_2, and the third inequality is by w_k= (A_k-A)(d_k+ x^s), k≥ N_1, and the definitions of D and ϵ_1. 
	
	
	It can be shown by induction that d_t ≤ D for any t ≥ N_1+N_2+1 in the same way, which completes the proof. 	

	
	






Prove d_t →0.
We will show that for any ϵ_3>0, there exists N_3, such that d_t_2≤ϵ_3 when t> N_3. The proof is very similar to the proof of Fact 2.

Since A_t→ A, when ϵ'_1= ϵ_3/(4c_6D), there exists N'_1, such that when t≥ N'_1, A_t-A≤ϵ'_1, where D is defined in Fact 2. 
Since A is stable, we have A^t  → 0, so when ϵ'_2= ϵ_3/(2D), there exists N'_2, such that when t>N'_2, A^t ≤ϵ'_2. 
Let N_3=N'_1+N'_2. When t>N_3,

    d_t+1   =  A^t-N'_1+1d_N'_1+ w_t +Aw_t-1+… + A^t-N'_1w_N'_1
       ≤A^t-N'_1+1 D +  w_t +Aw_t-1+… + A^t-N'_1w_N'_1_2
       ≤ϵ'_2 D + c_6 max_N'_1≤ k ≤ tw_k
       ≤ϵ'_2 D + 2c_6 ϵ'_1 D = (1/2+1/2)ϵ_3=ϵ_3

	where the second inequality is by Fact 1 and the definitions of ϵ_2, and the third inequality is by w_k= (A_k-A)(d_k+ x^s), k≥ N_1, max_k≥ 0(d_k, x^s)≤ D, and the definition of ϵ'_1. 

This completes the proof of the claim. 








 §.§ Proof of Lemma <ref>. 

Let x_t^*, u_t^* denote the optimal state and the optimal control input at time t respectively. By Theorem <ref>, the optimal controller is u_t^*=-K_t x_t^* +K_t^αα_t. For ease of notation, define

    D_t  A-BK_t.

Then, the dynamical system of x_t^* can be represented as 

    x_t+1^* = D_t x_t^* + BK_t^αα_t+1.


Proof outline:  We will prove x_t^* is bounded by three steps: 1) show that system x_t+1=D_t x_t is exponentially stable, 2) show that BK_t^αα_t+1 is bounded, 3) show x_t^* is bounded by the fact that exponentially stable systems are bounded-input-bounded-output stable.

Step 1: show x_t+1=D_t x_t is exponentially stable by a Lyapunov function.

Consider dynamical system x_t+1=D_t x_t.
	Define the state transition matrix: 
    Φ(t,t_0)=D_t-1⋯ D_t_0

	for t≥ t_0, and Φ(t,t)=I. 
	For any N, any 0≤ t_0 ≤ N t_0≤ t ≤ N, any Q_t ∈, R_t ∈, Q_N ∈, and for any x_t_0, the system is exponentially stable, i.e.
	
    Φ(t,t_0)    ≤ c_7 c_2^t-t_0
 
	where c_7=√(υ_max(P̅)/υ_min(P)), c_2=√(1-μ_f/υ_max(P̅))∈ [0,1).


We prove the exponential stability by constructing a Lyapunov function: L(t, x_t)=x_t^⊤ P_t x_t for t≥ 0. 

Claim:  For any x_t, the Lyapunov function satisfies

    υ_min(P)x_t_2^2 ≤ L(t,x_t)≤υ_max(P̅)x_t^2,   
    	 L(t+1, D_t x_t)-L(t,x_t)≤ -μ_f x_t_2^2.

	 where υ_max(·) and υ_min(·) denote the maximum and minimum eigenvalues of a matrix  respectively.
	 
Proof of Claim:
	By Lemma <ref>, P_t ∈, so
	υ_min(P) I_n ≤P≤ P_t ≤P̅≤υ_max(P̅)I_n.
	Thus, for any x_t, 
	
    υ_min(P)x_t_2^2 ≤ L(t,x_t)=x_t^⊤  P_t x_t≤υ_max(P̅)x_t^2

	
Besides,
	
    L(t+1, D_t x_t)-L(t,x_t)   = x_t^⊤ D_t^⊤ P_t+1D_t x_t -x_t^⊤  P_t x_t
        = x_t^⊤ (D_t^⊤ P_t+1D_t-P_t)x_t
       = x_t^⊤ (-Q_t-K_t^⊤ R_tK_t)x_t 
       ≤ -x_t^⊤Qx_t 
    	= -μ_f x_t^2

	where the third equality is by Theorem <ref>, the first inequality and the last equality are by Q_t+K_t^⊤ R_tK_t≥ Q_t ≥Q=μ_f I_n.
	By the claim above, 
	
    L(t+1, x_t+1)-L(t,x_t)≤ -μ_f x_t_2^2 ≤ -μ_f/υ_max(P̅)L(t,x_t)

	Thus, L(t+1, x_t+1) ≤ c_2^2 L(t,x_t) where c_2= √(1-μ_f/υ_max(P̅)). Here, c_2 is well-defined because 0≤μ_f I_n ≤Q̅≤P̅≤υ_max(P̅).
	
	
	
	For any t_t_0 and any x_t_0, it is easy to verify that the state x_t satisfies x_t=Φ(t, t_t_0)x_t_0. Therefore, 
	
    Φ(t, t_0)   = max_x_t_0≠0x_t/x_t_0
       ≤max_x_t_0≠0√(υ_max(P̅)/υ_min(P)L(t,x_t)/L(t_0, x_t_0))
       ≤√(υ_max(P̅)/υ_min(P)) c_2^t-t_0

	

Step 2: show that BK_t^αα_t+1 is bounded. We will show that 

    α_t ≤c_7/1-c_2υ_max(P̅)θ̅α̅,      BK_t^αα_t_2 ≤B_2^2 α̅/μ_g.
 

By Theorem <ref>, α_t satisfies the dynamical system α_t= D_t^⊤α_t+1+ Q_tθ_t, with initial condition α_N=Q_Nθ_N. As a result, we can write α_t in terms of θ_s and the transition matrix Φ(t, s) as follows

    α_t    =Q_tθ_t + D_t^⊤  Q_t+1θ_t+1+ …+ D_t^⊤… D_N-1^⊤  Q_Nθ_N
       =Φ(t,t)^⊤ Q_tθ_t+ Φ(t+1,t)^⊤  Q_t+1θ_t+1+ …+ Φ(N,t)^⊤  Q_Nθ_N

Then, the bound of α_t can be derived as follows.
	
    α_t   ≤Φ(t,t)^⊤Q_tθ_t+ …+  Φ(N,t)^⊤ Q_Nθ_N
        = Φ(t,t)Q_tθ_t+ …+  Φ(N,t) Q_Nθ_N
       ≤ c_7 υ_max(P̅) θ̅+ … +  c_7c_2^N-tυ_max(P̅)θ̅
       ≤ c_7 υ_max(P̅)θ̅1/1-c_2=α̅

	where the second inequality is by Lemma <ref>, Q_t ≤Q̅≤P̅ and Q_N ∈. 
	
	Consequently, 
		
    BK_t^αα_t    = B(R_t + B^⊤  P_t+1B)^-1 B^⊤α_t 
       ≤B^2 (R_t + B^⊤  P_t+1B)^-1α_t
       ≤B^2 α̅/μ_g

	




Step 3: bound x_t^*.

    x_t^*    = Φ(t,t)BK_t-1^αα_t+Φ(t,t-1)BK_t-2^αα_t-1+…Φ(t,1)BK_0^αα_1
       ≤ c_7 B^2 α̅/μ_g(1+c_2 + c_2^2 +…)= c_7 1/1-c_2B^2 α̅/μ_gx̅




















 §.§ Proof of Lemma <ref>

By Theorem <ref> and (<ref>), β _k = P_k^-1α_k≤1/υ_min(P)α̅ for any k, any N and any Q_t ∈, R_t ∈, Q_N ∈. When Q_t=Q ∈, R_t=R= for all t, β_k→β^e as k → -∞ by the proof (Part i)) of Proposition <ref>. Thus, β^e ≤1/υ_min(P)α̅. 
Define β̅=max(θ̅, 1/υ_min(P)α̅). This completes the proof.

 



§ SIMULATION DESCRIPTIONS
 



 §.§ LQT


The experiment settings are as follows. Let A=[0,1; -1/, 5/6 ],B=[0;1], N=30. Consider diagonal Q_t, R_t with diagonal entries i.i.d.  from Unif[1,2]. Let θ_t i.i.d. from Unif[-10,10]. 
 The stepsizes of RHGD and RHTM are based on the conditions in Theorem <ref>.  The stepsizes of RHAG can be viewed as RHTM with δ_c=1/l_c, δ_y=δ_ω=√(ζ)-1/√(ζ)+1 and δ_z=0. 




 §.§ Robotics tracking

Consider the following discrete-time counterpart of the kinematic model 


    x_t+1    = x_t + Δ t ·cosδ_t · v_t 
    y_t+1    = y_t+Δ t·sinδ_t · v_t
    δ_t+1    =δ_t + Δ t· w_t


Thus we have 


    δ_t     = arctan( y_t+1-y_t/x_t+1-x_t )
    
        v_t     = 1/Δ t·√((x_t+1-x_t)^2 + (y_t+1-y_t)^2)
    
        w_t    = δ_t+1 - δ_t/Δ t = 1/Δ t·[   arctan( y_t+2-y_t+1/x_t+2-x_t+1 ) - arctan( y_t+1-y_t/x_t+1-x_t )                 ]



So that (δ_t,v_t,w_t) can be expressed by the state variables (x_t,y_t).

In the simulation, the given reference trajectory is 


    x^r_t     = 16sin^3(t-6)
    
        y^r_t    = 13cos(t)-5 cos(2t-12)-2cos(3t-18)-cos(4t-24)


 As for the objective function, we set the cost coefficients as

    c_t = 
        0,     t = 0
    
        1,    otherwise    
            c_t^v = 
        0,     t = N
    
        15Δ t^2,    otherwise         c_t^w = 
        0,     t = N
    
        15Δ t^2,    otherwise

The discrete-time resolution for online control is 0.025 second, i.e., Δ t =0.025s. When implementing each control decision, a much smaller time resolution of 0.001s is used to simulate the real motion dynamics of the robot.


























