acker(): "Error using sym/poly". Why?

I'm trying to find the eigenvalues of a matrix H manually, where H represents some system dynamics.

(The system is defined by H, so eigs are solutions to det(H-lamba*I) = 0)

Once i find those eigs, i input them to acker() for pole placement to find gains.

(1) Why am i getting an error when inputting what i think are a set of normal numbers to acker()? See code below.

Note: I compare this to the normal eig(H) results to see if it's working, so i can see the normal method works with acker().

(2) Why is vpa() needed, ie why is solve() returning functions of z? I thought solve() would return numbers, since it's just trying to find solutions to det(H-lamba*I) = 0), which evaluates to:

lam^6 - 15999360000*lam^4 + 16000160000000000*lam^2 - 160000000000000000000000 == 0

and that seems like returned solutions should be just normal complex numbers.

A = [0 1; 
     -4e5 -400];
B = [0; 4e5];
C = [1 0];
Aaug = [A, zeros(2,1); -C, 0];
Baug = [B; 0];

H = 1.0e+12 * [0   0.000000000001000                   0                   0                   0                   0;
  -0.000000400000000  -0.000000000400000                   0                   0  -0.160000000000000                   0;
  -0.000000000001000                   0                   0                   0                   0                   0;
  -0.000000100000000                   0                   0                   0   0.000000400000000   0.000000000001000;
                   0  -0.000000000000100                   0  -0.000000000001000   0.000000000400000                   0;
                   0                   0  -1.000000000000000                   0                   0                   0];

% ~~~~~~~ Built-in eigs ~~~~~~~
% These are a set of numbers
eigs1 = eig(H)

eigs1 = 
   1.0e+05 *

  -1.2648 + 0.0000i
   1.2648 + 0.0000i
  -0.0135 + 0.0115i
  -0.0135 - 0.0115i
   0.0135 + 0.0115i
   0.0135 - 0.0115i

whichNegEig = find(real(eigs1)<0);
K = acker(Aaug, Baug, eigs1(whichNegEig))

K = 1×3

1.0e+06 * 0.0009 0.0000 -1.0000

% ~~~~~~~ Manually find eigs ~~~~~~~
syms lam;
sysDet = det(H - lam*eye(size(H)));

polySys = solve(sysDet == 0, lam)

polySys = 

eigs2= vpa(polySys)

eigs2 = 

% These SHOULD be a set of numbers. Are they not?
whichNegEig2 = find(real(eigs2)<0);
% They SEEM the same as whichNegEig, so seem correct.
eigs2(whichNegEig2)

ans = 

% Why is this giving this error:
% "Error using sym/poly
% SYM/POLY has been removed. Use CHARPOLY instead."

K2 = acker(Aaug, Baug, eigs2(whichNegEig2))

Error using sym/poly
SYM/POLY has been removed. Use CHARPOLY instead.

Error in acker (line 41)
k = ctrb(a,b)\polyvalm(real(poly(p)),a);

[John]

Huh -- per this post, it looks like double() is needed, ie

acker(Aaug, Baug, double(eigs2(whichNegEig2)))

With double(), the error goes away. Why is that needed?

Also, question (2) still stands.

[Star Strider]

The Symbolic Math Toolbox (SMT) calculations are always all symbolic, not numeric in the sense that they can be used outside the SMT . The vpa function converts the fractional representation (including ‘root’) into their numeric equaivalents, however they remain symbolic variables. So the vpa call is not absolutely necessary, however double is, because it converts the symbolic results to numeric results that can be used elsewhere in MATLAB calculations..

.

[John]

Thanks @Star Strider. "it converts the symbolic results to numeric results" How can i tell they are symbolic results?

IeI tried a few "isxyz()" commands -- how would i identify it as a symbolic type?

I also tried

sympref('FloatingPointOutput',true);

before solve(), but no change. I thought that would give floating point outputs, if possible.

[Star Strider]

My pleasure!

If it is not obvious by the context, one option is to use the whos function.

Example —

x = 1

x = 1

whos x

  Name      Size            Bytes  Class     Attributes

  x         1x1                 8  double              

syms x
y = x

y = 

x

whos y

  Name      Size            Bytes  Class    Attributes

  y         1x1                 8  sym                

The sympref call:

sympref('FloatingPointOutput',true);

produces floating-point representation with four digits after the decimal, instead of producing fractional constants. It does not automatically invoke vpa. See Display Symbolic Results in Floating-Point Format for a full explanation.

.

Answer 1

[Paul]

Hi John,

A = [0 1; 
     -4e5 -400];
B = [0; 4e5];
C = [1 0];
Aaug = [A, zeros(2,1); -C, 0];
Baug = [B; 0];
H = 1.0e+12 * [0   0.000000000001000                   0                   0                   0                   0;
  -0.000000400000000  -0.000000000400000                   0                   0  -0.160000000000000                   0;
  -0.000000000001000                   0                   0                   0                   0                   0;
  -0.000000100000000                   0                   0                   0   0.000000400000000   0.000000000001000;
                   0  -0.000000000000100                   0  -0.000000000001000   0.000000000400000                   0;
                   0                   0  -1.000000000000000                   0                   0                   0];
syms lam;
sysDet = det(H - lam*eye(size(H)));

Here, solve trying to solve for closed form expressions of a sixth order polynomial, which can't be done in general. In some cases, solve will revert to vpasolve, but I guess not in this case.

polySys = solve(sysDet == 0, lam)

polySys = 

However, here sysDet has a form that solve can (arguaby) handle a bit better if we ask it to

polySys = solve(sysDet == 0, lam,'MaxDegree',3)

polySys = 

At this point, polySys is of class sym

class(polySys)

ans = 'sym'

Now, we can convert polySys to "numbers"

eigs2= vpa(polySys)

eigs2 = 

but eigs2 is still of class sym

class(eigs2)

ans = 'sym'

Calling acker with sym inputs just isn't allowed, hence the need to convert to double. Some functions outside the Symbolic Math Toolbox that ordinarily expect numeric inputs can deal with sym inputs, but acker isn't one of those.

Having said that, acker is pretty simple and I think you could implement it pretty easily with sym variables, at least the basic part that computes the gains.

[Walter Roberson]

solve() returns in terms of roots() for polynomials. solve() only invokes vpasolve() for nonlinear systems.

Note: vpasolve() returns all numeric roots of polynomials that are the same degree as the number of variables asked to solve for. This is an exception to the more general rule that vpasolve() returns only a single numeric value.

[John]

Ah, good to know, thanks @Walter Roberson. That's helpful.

[John]

Thanks @Paul -- that helps me understand better what options are available.

'MaxDegree',3

does indeed seem to convert directly to values (syms), and not root equations. But...why does that work?

Presumably, the order of solve()'s solution is already known to Matlab, and since it's <5 (per the documentation) it should be calling an explicit representation solution.

[Paul]

I don't know why the default for MaxDegree is 2. My wild guess is that the expressions for the roots of higher order polynomials are pretty ugly and not particularly useful or insightful in general, particularly when the coefficients of the polynomial themselves include sympolicy variables. But that's just a guess.

Also, if you don't mind me asking, why use symbolic here at all? All that it seems to offer in this particular example is the extended precision in the computation of eigs2.

[Walter Roberson]

The explicit form of quadratic roots is easy to understand. The explicit form of cubic roots is somewhat messy and I have met very few people who understand it fairly directly. I have never encountered anyone who understood the explicit form of quartics directly: for example if you were to flip the sign of one of the terms would they even notice or understand the significance?

[John]

@Paul "Also, if you don't mind me asking, why use symbolic here at all? All that it seems to offer in this particular example is the extended precision in the computation of eigs2."

That's a great question: I don't actually need symbolic, but i didn't know of another straightforward method to manually compute eigenvalues for a 6th-0 equation, other than using poly solve via solve(). Do you have any suggestions?

@Walter Roberson 100% agreed -- but then why not return the number, vs an expression, and abstract away the underlying complexity?

[Paul]

I'm not following. eigs2 are the eigenvalues of H, which could just as well be computed with

eig(H)

albeit only in double precision.

[John]

Sorry : i meant a manual computation per the original post's motivation:

'I'm trying to find the eigenvalues of a matrix H manually", etc

Eig() is not exportable to C/C++, so i'm trying to find other methods... solve() is a step closer.

[Paul]

It isn't? According to its doc page, eig is eligible for C/C++ code generation. I don't do code generation, so don't know any more than that.

[Walter Roberson]

Nothing in the Symbolic Toolbox can be compiled with MATLAB Compiler or MATLAB Coder.

roots can be used to compute numeric roots of a polynomial that has been abstracted to its list of coefficients of descending powers.

A strategy that can be used is to use an interactive session to generate the symbolic eigenvalues down to the root() or rootOf() structure. Then use children() to extract the symbolic polynomial from inside the root() or rootOf() structure, and use sym2poly() to get the list of coefficients of each power, in symbolic form.

Now you can matlabFunction() that list, asking to write to a file; you might need to request Optimize false unless they have fixed the optimization bugs from the last several releases.

The result would be a .m file that you could invoke from a different program -- the program that needs to be compiled.

The symbolic form of eigenvalues is often pretty long, and most of the time it is faster and more accurate to use eig() -- which has been supported for code generation for a number of releases (earlier than R2018b)

Answer 2

[Walter Roberson]

acker() is a now-undocumented function in the Control Systems Toolbox, in a directory reserved for "obsolete" functions. In other words, it was withdrawn and should not be used in any new code.

I see people recommending place instead.

[John]

Thanks @Walter Roberson, will change.

[Paul]

Keep in mind that place has a constraint on the desired pole locations.

rng(100)
A = rand(2);B = rand(2,1);
place(A,B,[-2 , -2]);

Error using place
The "place" command cannot place poles with multiplicity greater than rank(B).

[John]

Yes, thanks. That's why i was hoping to use acker(). But, i guess this is now a limitation.