Symbolic Solutions That Populate a Vector

I am trying to approximate a solution to the Riccati equation. To do this, I have a final value and work backwards to obtain all previous values. The issue is that in this particular case the solution needs to be symbolic bexause it depends on a function of time that will be calculated later. The code I came up with is listed below. This works very well when there is no symbolic variable to deal with.

% P(tf)
P_11(length(t)) = F(1,1);
P_12(length(t)) = F(1,2);
P_22(length(t)) = F(2,2);

% Solve Riccati eqn backwards
syms z;
for n=length(t) : -1 : 2
    P_11_dot = P_12(n)^2 + 3*z*P_12(n) + 3*P_12(n)*z - 1;
    P_12_dot = 2*P_12(n) - P_11(n) + P_12(n)*P_22(n) + 3*P_22(n)*z;
    P_22_dot = 4*P_22(n) - 2*P_12(n) + P_22(n)^2;
    P_11(n-1) = P_11(n) - P_11_dot*step;
    P_12(n-1) = P_12(n) - P_12_dot*step;
    P_22(n-1) = P_22(n) - P_22_dot*step;
end

the error message i get is this;

The following error occurred converting from sym to double:

Unable to convert expression into double array.

Error in LTV_SYS (line 56)

P_11(n-1) = P_11(n) - P_11_dot*step;

The idea is to then calculate z using the code below. Here, I know the initial conditions and since my symbolic variable (z) is a state, I should be able to substitute in for it. The problem is I can't get to this point because of the above error.

% Initialize state vector
x(:,1) = x0;

% Solve for optimal state and control
for n=1 : length(t)-1
    z = x(1,n);
    A = [0 1; -3*z -2];
    
    P = [P_11(n) P_12(n); P_12(n) P_22(n)];
    K(n,:) = inv(R)*B'*P;
    U(n+1) = -K(n,:)*x(:,n);
    x_dot = A*x(:,n)+B*U(n+1);
    x(:,n+1) = x(:,n) + x_dot*step;
end

I appreciate any help that anyone can offer.

Answer 1

[Walter Roberson]

P_11(length(t)) = sym(F(1,1));
P_12(length(t)) = sym(F(1,2));
P_22(length(t)) = sym(F(2,2));

[Justin]

Hi Walter,

Thank you for your reply. I did as you suggested and initially everything seemed okay. However, after nearly 2 hours of time I stopped the code from executing. It looks to me like the coefficients may be taking up too much memory. After 10/1000 iterations I get the following warning;

Output truncated. Text exceeds maximum line length for Command Window display.

Essentially, I believe my general approach may not be ideal. It is far too computationally exhaustive. Can you (or anyone else) suggest a more efficient alternative?

Thanks,

Justin

[Steven Lord]

That message is not a warning. It indicates that the text you're trying to display is longer than 32,768 characters (I believe that's the current threshold.) The rule of thumb I've seen is that one page of text is about 3,000 characters so your expression is longer than 11 pages or so. Do you really need to read through an expression that is longer than 11 pages? [Actually, I did a computation that triggered this message and pasted the result into Microsoft Word. With normal margins and 11 point font it just about filled 9 pages.]

A couple more things I noted looking at your original code:

% Initialize state vector
x(:,1) = x0;

% Solve for optimal state and control
for n=1 : length(t)-1

Roughly how long do you expect t to be? Is length(t) more like ten, ten thousand, ten million, etc.?

    z = x(1,n);
    A = [0 1; -3*z -2];
    
    P = [P_11(n) P_12(n); P_12(n) P_22(n)];
    K(n,:) = inv(R)*B'*P;

Is that a typo (and you intended to invert P) or are you really inverting the (symbolic?) matrix R that doesn't change inside the loop at each iteration? If the latter, compute that symbolic inverse once before the loop and use it inside the loop. [Actually, if it's not a typo compute inv(R)*B' outside the loop and make this step just a matrix-matrix multiplication, since B also does not change inside the loop.]

    U(n+1) = -K(n,:)*x(:,n);
    x_dot = A*x(:,n)+B*U(n+1);
    x(:,n+1) = x(:,n) + x_dot*step;
end

[Justin]

Hi Steven,

Thanks for your reply.

In this case, length(t)=1001. There is no typo regarding inv(R). R is a weight matrix that is constant and does not vary with time. I could certainly calculate inv(r)*B' outside of the loop but the k variable is a gain and it is important for me to know its value at all time (t).

The problem I have is that P is a function time and is dependent on z which is also a function of time. I have final values for P and need to work my way backward to the beginning of time to solve for P(1). This is a typical approach when finding solutions to the differential Riccati equation. What makes this problem untypical is the dependence of P on z. I have initial values for z and can work forward in time once I know P since z(t + ) is dependent on P(t).

Any ideas on how to compute in a reasonable amount of time? My current approach is failing miserably.

[Walter Roberson]

At the moment, you are computing indefinite precision rational solutions. That is fine in theory if your initial F values are indefinite precision, but we can see from the error message that they are floating point numbers (double). Unless they are integers expressed as floating point, the inheirent uncertainty in the initial values is going to get magnified by the calculation, and unless you were intending to do an error analysis watching how the solution changes as you change the uncertainty, it becomes difficult to justify creating huge rational expressions.

You should probably consider using vpa() at appropriate locations to keep the length of the expressions under control. You can use digits() to control the precision you want to execute to.

For example I might convert

    K(n,:) = inv(R)*B'*P;
    U(n+1) = -K(n,:)*x(:,n);
    x_dot = A*x(:,n)+B*U(n+1);
    x(:,n+1) = x(:,n) + x_dot*step;

into

    this_K = inv(R)*B'*P;
    this_U = -this_K*x(:,n);
    x_dot = A*x(:,n)+B*this_U;
    K(n,:) = vpa(this_K);
    U(n+1) = vpa(this_U);
    x(:,n+1) = x(:,n) + vpa(x_dot*step);
    

This preserves intermediate precision.

[Justin]

The current F matrix is a diagonal matrix comprised of only integers (). I only need 3 of the elements of F. The elements of F used for P at time final are a 4 and two zeros.

Currently the second loop (where I solve for optimal control and states) is commented out. I cannot get out of the first loop where I "solve" (really approximate) the Riccati equation.

I am not going to do any error analysis, however, I do neet to know P_xy at all times. Since my attempt is an approximation, I agree that I do not need this degree of precision. You are correct, is not justified. Do you have any advice on how I can deal with it?

[Justin]

I tried using vpa in the P loop as shown below. It is now approximating the rational coefficients as decimals. This seems like a good start but I think it is still providing too much precision which is slowing down computation. Is there a way to round the results to, say 5, decimal places?

for n=length(t) : -1 : 2
    P_11_dot = vpa(P_12(n)^2 + 3*z*P_12(n) + 3*P_12(n)*z - 1);
    P_12_dot = vpa(2*P_12(n) - P_11(n) + P_12(n)*P_22(n) + 3*P_22(n)*z);
    P_22_dot = vpa(4*P_22(n) - 2*P_12(n) + P_22(n)^2);
    P_11(n-1) = vpa(P_11(n) - P_11_dot*step);
    P_12(n-1) = vpa(P_12(n) - P_12_dot*step);
    P_22(n-1) = vpa(P_22(n) - P_22_dot*step);
end

[Walter Roberson]

"You can use digits() to control the precision you want to execute to."