nonlinear variable state space observer pole placement

Hi,

i want to place some observer poles in dependency of an variable. The variable v is speed and it is changing with the time.

Example:

1. option

syms v

A = [-1 2; v -1]

c = [0 1]

p = [-2; -3];

so,

L = place(A',c',p).'; %doesnt work

-----------

2.option

syms v lambda l1 l2

L =[l1;l2];

solve((det(lambda * [1 0; 0 1] - (A - L*c)) == (lambda+2)*(lambda+3)),[l1 l2]); %matlab doesn't eliminate lambda

l1 = ?

l2 = ?

do i need to calculate the det matrix by myself or can matlab help me in a faster way?

Thanks and Regards

Bernd

Answer 1

[Paul]

Hi Bernd,

Matlab assumes by default that all sym variables are complex, and this assumption can sometimes lead to unexpected results.

Your code works if we explicitly assume that v and the estimator gains are real (my experience is to include all relevant assumptions)

syms v l1 l2 real
syms lambda
A = [-1 2; v -1];
c = [0  1];
p = [-2; -3];
L = [l1;l2];
sol = solve((det(lambda * [1 0; 0 1] - (A - L*c)) == (lambda+2)*(lambda+3)),[l1 l2])         %matlab doesn't eliminate lambda

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

sol = struct with fields:
    l1: (2*v + 2)/v
    l2: 3

If v == 0, then A,C are not an observable pair, so we should exclude v == 0 from the analysis

assumeAlso(v ~= 0);
sol = solve((det(lambda * [1 0; 0 1] - (A - L*c)) == (lambda+2)*(lambda+3)),[l1 l2])

sol = struct with fields:
    l1: (2*v + 2)/v
    l2: 3

Or more compactly without hardcoding anything

sol = solve(det(lambda*eye(2)- (A - L*c)) == poly2sym(poly(p),lambda),L)

sol = struct with fields:
    l1: (2*v + 2)/v
    l2: 3

Or, using Ackerman's formula

alpha(lambda) = poly2sym(poly(p),lambda)

alpha(lambda) = 

O = [c ; c*A];
(A^2 + 5*A + 6*eye(2))*inv(O)*[0;1]

ans = 

[Bernd Pfeifer]

Hi Paul,

thanks for the answer, that already helped a lot.

It worked with the dummy matrices, but with the actual larger matrices I still get a problem:

Here are the actual matrices and code:

syms v l11 l12 l13 l14 l21 l22 l23 l24 real

syms lambda

A = [-4, v, 4, 0;

-v,-4, 0, 4;

0, 0, 0, 0;

0, 0, 0, 0];

C = [10 0 -10 0

0 10 0 -10];% system is obeservable with rank([C ; C*A; C*A*A; C*A*A*A]) = 4

L1 = [l11,l12;l13,l14];

L2 = [l21,l22;l23,l24];

L = [L1;L2];

assumeAlso(v ~= 0);

sol = solve((det(lambda * eye(4) - (A - L*C)) == poly2sym(poly(p),lambda)),[l11 l12 l13 l14 l21 l22 l23 l24])

It still has the lambdas in the solution. Do you have a solution for that? :)

Thanks a lot!

Regards

Bernd

[Bernd Pfeifer]

ah sry, i forgot the poles:

p = [-5; -5; -10; -10].*1000;

[Paul]

It still has lambdas in the solution because we now have an undertermined system, i.e., four equations with eight unknowns. So you'll need to either specify four of the observer gains a priori or include four other equations (or a combination of both) to uniquely determine the gains. This sounds easier than it is.

For state feedback controller design, there are methods for the multivariable case for partial eigenstructure assigment. By partial here I mean that there is never enough degrees of freedom to fully specify the closed-loop eigenvectors and the closed-loop eigenvalues, but we can specify the eigenvalues and then get eigenvectors that are close, in some sense, to desired eigenvectors. I see no reason why a similar technique couldn't work for observer design. There may be other methods as for muti-output observer design, such as specifying the design parameters for a Kalman filter and then taking the steady-state gains, for example. Or the algorithm in the Reference [1] of place. Of course, implementing any of these techniques symbolically may be a chore.

[Paul]

Here is one approach

syms v real
A = [ -4, v, 4, 0;
      -v,-4, 0, 4;
       0, 0, 0, 0;
       0, 0, 0, 0];      

C = [10    0  -10    0
      0   10    0  -10];% system is obeservable with rank([C ; C*A; C*A*A; C*A*A*A]) = 4

This example works with repeated poles in p, but I'm not sure it will in general, so I'll change p so that the estimator poles are distinct.

%p = [-5; -5; -10; -10].*1000;
p = [-5; -6; -10; -11].*1000;

rng(100)
for ii = 1:4
    Q(:,ii) = rand(2,1);
    V(:,ii) = inv(A.' - p(ii)*eye(4)) * C.' * Q(:,ii);
end

L = (Q*inv(V)).';

sort([eig(A - L*C) p])

ans = 

Choosing estimator poles so far out into left half plane could be problematic, both from a computational perspective and from using such a high bandwidth design if using this estimator for closed-loop control.

[Bernd Pfeifer]

Okay true i understood the problem.

But i didn't understood the outcome of your example approach. Which technique was this and how should i implement it so solve the problem?

I am now trying to find a solution in the Reference [1] of place, where "the algorithm uses the extra degrees of freedom to find a solution that minimizes the sensitivity of the closed-loop poles to perturbations in A or B."

[Paul]

I'm not sure the technique has a name. It's basically using the relationship between the eigenvalues and left eigenvectors of (A - L*C).

I don't understand the question "how should I implement it"? Doesn't the code show how to implement it? The entries in L are functions of v. For any value of v the observer gain matrix L combined with A, both evaluated at v, and combined with C, result in estimator poles of A - L*C at the specified locations.

[Bernd Pfeifer]

okay the implementation is clear.

But did i understood it right? You calculated the eigenvectors with (A-p*E)v = C' *Q?

But where do C' *Q come from? Do you give the equation some random values, because the system of equations is undertermined?

[Paul]

Each V_i is a left eigevector of (A-LC) that corresponds to p_i.

Then

So the idea is to specify the 2x1 Qi, i = 1-4, to compute the 4x1 Vi, which are the left eigenvectors of (A - LC).

I gave the Q_i random values because the problem gave no specification on what the Vi should be. It may be possible to specify a desired set, call it V_d, then for each V_d_i find the Q_i that yields a V_i that is in some sense close to V_d_i, then proceed as above with the full set of Q_i and V_i to solve for L.

Try searching on state feedback for eigenstructure assignment. Might get most of the hits to solve for the control problem, but any of those should be adaptible to the estimation problem.

[Bernd Pfeifer]

Perfect, i understood the most now and the solution works!

Thanks a lot!

Best Regards

Bernd

[Bernd Pfeifer]

One question left.

Is there a reason why you chose the left eigenvectors for the assignement and not the right eigenvectors?

Regards

Bernd

[Paul]

Assuming you mean the right eigenvectors of A - LC, I didn't see how that could work, because we need L (or transpose(L) ) to multiply V, but with a right eigvector we get a term LCV.

Answer 2

[Sam Chak]

Hi @Bernd Pfeifer

Long time I didn't solve math puzzles like this one. Generally, you should be able find the analytical solution for L that satisfies the condition of eigenvalues of that gives and .

The following is an attempt to solve the problem heuristically for . You can attempt for .

v = 1:10;
c = [0 1];
p = [-2; -3];
for j = 1:length(v)
    A = [-1 2; v(j) -1];
    L = place(A', c', p)   % always return L = [m 3]
    m(j) = L(:, 1);
end

L = 1×2

4.0000 3.0000

L = 1×2

3.0000 3.0000

L = 1×2

2.6667 3.0000

L = 1×2

2.5000 3.0000

L = 1×2

2.4000 3.0000

L = 1×2

2.3333 3.0000

L = 1×2

2.2857 3.0000

L = 1×2

2.2500 3.0000

L = 1×2

2.2222 3.0000

L = 1×2

2.2000 3.0000

plot(v, m, 'p'), xlabel('v'), ylabel('L_1')

The pattern can be intuitively recognived as:

vv = linspace(1, 10, 901);
L1 = (4 + 2*(vv - 1))./vv;
plot(vv, L1), grid on, xlabel('v'), ylabel('L_1')

[Bernd Pfeifer]

Hi Sam,

thanks for the answer, but i want to implement the observer in simulink so i need one equation.

Paul already helped a lot and i think we get the problem in his way solved :)

Regards

Bernd

[Sam Chak]

@Bernd Pfeifer, don't mention it. @Paul is a helpful and knowledgeable person in control design.

For implementation in Simulink, both equations for are exactly the same (for your original 2nd-order system):

My equation came from my recognition of the numerical pattern as something related to the arithmetic sequence (without employing the curve-fitting tool), while Paul's equation is the solution as a direct result from his analytical mind and the computational power of MATLAB.