1

Is it as straightforward to add time delay compensation to LQR that uses full-state feedback, as it is with PID that uses error-based feedback?

PID can use eg Smith predictor, but when I search for time delay compensation with LQR, I basically only see research papers; is there any good standard technique for that?

Maybe I don't know what terms to search for, but "LQR with time delay compensation" yields minimal actionable results.

J B
  • 35
  • For LQR it is usually assumed that the full state is known. So in your case only a time delayed full state is known? And are you considering continuous or discrete time? – Kwin van der Veen May 06 '24 at 05:58
  • Correct, assume the full delayed state is known, as well as the time delay. I'm considering the discrete case, but continuous is fine too for now -- essentially, i'm looking to understand the basic principle, since the lqr architecture is different than an error-based controller. – J B May 06 '24 at 06:36

2 Answers2

3

I assume here that you consider a system with input delay of the form

$$\dot{x}(t)=Ax(t)+Bu(t-h)$$

where $h>0$ is the delay. A state-predictor is given in this case by

$$ x_p(t):=e^{\tilde A\tilde h}x(t)+\int_0^{\tilde h}e^{\tilde As}\tilde Bu(t-s)ds $$ where $\tilde A,\tilde B$ and $\tilde h$ are the known values for $A,B$, and $h$. When there are no uncertainties, i.e. $\tilde A=A$, $\tilde B=B$, and $\tilde h=h$, we have that $x_p(t)=x(t+h)$.

Using then perfect knowledge of the system, one can consider the control law $u(t)=Kx_p(t)$ to get the closed-loop system $$ \dot{x}(t)=(A+BK)x(t), $$ for which stability is easy to enforce assuming that $(A,B)$ is stabilizable. This procedure is called the Finite Spectrum Assignment (FSA). Now, you can use your favorite design method for finding the state-feedback matrix. The difficult part now is the implementation of the predictor. This is a well-studied problem as the naive approximations may lead to an unstable behavior.

KBS
  • 7,903
  • Thank you, I appreciate the details and help. This makes sense, but here's the step that i don't follow: "Now, you can use your favorite design method for finding the state-feedback matrix". How would i find K (state feedback matrix) when xp is of a form that requires numeric analysis? Commands like matlab or python lqr() require linear matrix inputs, so perhaps this isn't doable without manual solving of A.R.E, eg via Schur method. – J B May 21 '24 at 19:22
  • (2) I tried using x^ = exp(T(A-BK)) y^ for LQR using u = -Kexp(T(A-BK)), but that did not stabilize the system. What might i be doing wrong? – J B May 21 '24 at 19:44
  • You find $K$ such that $A+BK$ is stable and then you need to implement the predictor and the state-feedback. I am not sure how you came up with the exponential, but this is not how it can be done. You need to consider the predictor. – KBS May 22 '24 at 10:32
  • Thanks.Yes, it's clear that A-BK needs to be stable. The exponential form was what Kwin mentioned in their reply, referencing your reply. I thought the exponential WAS the predictor: isn't x * exp(T(A-BK)) the future state of x after time T, given the state evolution exp(T (A-BK)) ? Since state evolve per exp(DynamicsMatrix * t). If this isn't the way, how would i practically implement the predictor? – J B May 22 '24 at 23:16
  • The predictor is that integral expression I gave. There are ways to practically implement that expression using "good" discretization approaches. – KBS May 23 '24 at 04:10
  • I see. Thanks, makes sense to use the full integral expression (I must have misunderstood Kwin's comment). In any case, would you be able to point me to a practical "good" implementation? I'm not sure how to find one. By the way: interestingly, book-pg.279 (eq B.9) of this link ALSO shows the exp(T*A-BK) as a forward state propagation for latency compensation; though when I tried it, it didn't seem to work. https://file.tavsys.net/control/controls-engineering-in-frc.pdf. – J B May 23 '24 at 06:50
  • Check the book by Zhong, "Robust Control of Time-Delay Systems" for discussions about the implementation of the state predictor. – KBS May 23 '24 at 07:45
  • Great, thank you. – J B May 23 '24 at 08:06
1

For the Linear Quadratic Regulator problem the state is considered known and only really considers the optimal input that minimizes a specific type of cost function. Adding delay into the the mix turns the problem into the combination of state estimation and full state feedback. Due to the certainty equivalence principle for LTI systems, one can separate the state estimation from the optimal full state feedback.

The system for continuous time can be formulated as

\begin{align} \dot{x}(t) &= A\,x(t)+B\,u(t), \tag{1a} \\ y(t) &= C\,x(t-\tau), \tag{1b} \end{align}

with $\tau>0$ and $C=I$, but for generality I will keep using $C$, while assuming $(A,C)$ is observable. The full state feedback from LQR will be off the form $u(t) = -K\,\hat{x}(t)$, with the $\hat{x}(t)$ the state obtained from the state estimator. For example for the state estimation one can use a Luenberger observer,

\begin{align} \dot{\hat{x}}(t) &= A\,\hat{x}(t) + B\,u(t) - L(t)\,(\hat{y}(t) - y(t)), \tag{2a} \\ \hat{y}(t) &= C\,\hat{x}(t-\tau). \tag{2b} \end{align}

The separation due to the certainty equivalence principle can be obtained using $\tilde{x}(t) = x(t) - \hat{x}(t)$ giving

\begin{align} \dot{x}(t) &= (A - B\,K)\,x(t) + B\,K\,\tilde{x}(t), \tag{3a} \\ \dot{\tilde{x}}(t) &= A\,\tilde{x}(t) - L\,C\,\tilde{x}(t-\tau), \tag{3b} \end{align}

which is stable as long as $A-B\,K$ is Hurwitz and the autonomous dynamics of $\tilde{x}(t)$ is stable. Assuming that there is a $M$ such that $\tilde{x}(t-\tau) = M\,\tilde{x}(t)$ yields

$$ \dot{\tilde{x}}(t) = (A-L\,C\,M)\,\tilde{x}(t). \tag{4} $$

Solving this autonomous linear differential equation gives

$$ \tilde{x}(t-\tau) = \underbrace{e^{-\tau\,(A-L\,C\,M)}}_M\,\tilde{x}(t), \tag{5} $$

confirming the assumption that there is a $M$ such that $\tilde{x}(t-\tau) = M\,\tilde{x}(t)$.

Since $M$ is a solution to matrix exponential, means that for finite $\tau$ is holds that it is invertible. Furthermore, in your case we have $C=I$, which is also invertible. Therefore, in the general case, when $C$ is invertible, then $L$ can be chosen such that $A-L\,C\,M$ can be set equal to any Hurwitz matrix $X$ we want, using

$$ L = (A - X)\,M^{-1} C^{-1}. \tag{6} $$

This also means that $M$ and $M^{-1}$ can be solved for using

\begin{align} M &= e^{-\tau\,X}, \tag{7a} \\ M^{-1} &= e^{\tau\,X}. \tag{7b} \end{align}

In conclusion for the continuous time case one can stabilize $(1)$ using $u(t) = -K\,\hat{x}(t)$ and $(2)$, to solve for $\hat{x}(t)$, with $K$ for example obtained from LQR and $L$ obtained by choosing a Hurwitz matrix $X$ and substitute it in $(7)$ for $M$ and substitute $X$ and $M$ in $(6)$ to solve for $L$.

It can be noted that for non-invertible $C$ this procedure becomes a little bit more complicated, because solving for $L$ and $M$ becomes more intertwined.

Similar to the answer by KBS for the case of $C=I$ it is also possible to solve the delayed state forward in time, just as in $(5)$. In that case no state estimation is needed and instead $\hat{x}(t)=x(t)$ can be obtained, given $y(t)=x(t-\tau)$, using

$$ \hat{x}(t) = e^{\tau\,(A-B\,K)} y(t). \tag{8} $$

Kwin van der Veen
  • 12,770
  • Thank you, this is clearly laid out and helpful conceptually. To make sure I understand: you're saying solve via LQR using standard A,B,C (assuming standard no delay sys). But this is stabilized by the state estimator including the delay to give the appropriate x,y feedback, with L being found by a choice of X. The concept is that the system will be stabilized with standard gains as long as the state feedback is correctly estimated including delay. Is this right? If so, how is X chosen? That would be the one new input design matrix. – J B May 21 '24 at 19:31
  • @JB common rule of thumb is that the state estimator has about an order of magnitude faster exponential decay as the full state feedback. So if $\lambda$ is the real part of the eigenvalue of $A-B,K$ which has the most negative real part, one could choose $X=10,\lambda,I$. – Kwin van der Veen May 21 '24 at 20:26
  • I see -- thanks, that makes sense. Using the above i get L<0, though: since X is by design > 10lambdaI, isn't (A-X) < 0? I used an example A corresponding to an LQR problem, and used the closed-loop eigs to find lambda, and then X. That X is much larger than A. Am i misunderstanding (A-X) M^-1 C^-1 ? – J B May 22 '24 at 00:48
  • Ah -- you'd said to take the real part of the (negative) eig, not the magnitude. So, that makes X negative, and A-X is positive. Two questions as i'm working through your approach: (a) In the case of a system with 2 states and an A that's 2x2, how is C^-1 possible to use alongside an M^-1 that's 2x2? My system has C = [1 0], ie 1x2, so y = Cx = 1x1. (b) My M^-1 is ~ 1e-200, ie --> 0, and L --> 0. I use matlab's expm(tauX), with X = [-240 0; 0 -240]. But aren't matrix exponentials with negative matrix elems from tauX, bounded by +1? If so, L would also be small in what you showed. – J B May 22 '24 at 07:18
  • Typo above. I meant: "tau*X = [-240 0; 0 -240]" – J B May 22 '24 at 08:15
  • @JB Note that $\lambda$ is negative. And note that $\lambda$ is obtained from the eigenvalues of $A-B,K$, so nothing can be said about the eigenvalues of $A-L$ in general. – Kwin van der Veen May 24 '24 at 10:53
  • That makes sense. (a) How is (A-X) M^-1 C^-1 possible to compute, in the case that C is (1 0)? The dimensions don't seem to match. (b) Also, when i compute M^-1, it's ~1e-200; that would make L --> 0 in equation (6). Isn't M^-1 always forced to a small number, like what i was seeing, since X has negative elements per your comment above? – J B May 25 '24 at 04:58
  • @JB as mentioned in my answer I assume $C=I$, but this also works when $C$ is full rank. So for your $C$ the equations 6 and 7 don't work and you would have to solve $L$ and $M$ at the same time, which is more complicated. – Kwin van der Veen May 25 '24 at 08:09
  • Ah, quite so, thanks -- I'd read that several times and then not put it into context. Thanks :) I'm not sure how to solve both at the same time, but will ponder. Aside from that, I'm curious if you have any thoughts about question (b) from my last post: isn't M^-1 always be forced to a small number, given how X is found via M^-1 = e^(TX) since X has negative elements? – J B May 26 '24 at 01:32
  • @JB Correct, but do note that this does assume again that $C=I$ or at least invertible. In general it does not hold that all elements of $M^{-1}$ are bound by +1, just its eigenvalues. From a mathematics point of view $M$, $M^{-1}$ and $L$ can be calculated, but when doing so using numerical calculations on a computer one might run into numerical accuracy issues. – Kwin van der Veen May 26 '24 at 09:45
  • Thanks Kwin. For sake of my understanding $M$, assume C = I (though in my case it's not). Makes sense that "...it does not hold that all elements of $M^{−1}$ are bound by +1, just its eigenvalues". But, the eigenvalues of $M$ help determine $L$. If X=10λI, and λ<0, then eigs of T*X<0. Wouldn't that cause eigs of $M^{−1}$ to be very small numbers, so $L \rightarrow 0$? (I.e. observer doesn't converge). My example is λ = -3e6 rad/s, and T = 10 us: this causes $M^{−1}$ to have eigenvalues ~1e-138. Then L ~= 0. How can L be a reasonable magnitude if $L=(A−X) M^{−1} C^{−1}$, when $M^{−1}$ << 1? – J B May 26 '24 at 22:05
  • @JB because $M$ becomes very large, thus the order of magnitude of the elements of $L,C,M$ remain well defined. But the individual matrices do become less well conditioned. This is an inherent limitation of time delay, which is similar to right-half-plane zeros as demonstrated here. – Kwin van der Veen May 26 '24 at 23:54