Annoying VPA Bug that Generates Wrong Data Phase or Wrong Imaginary Value in a Complex Number

Run the following script and compare the numeric data and plots to see the VPA bug. This is tested with Matlab 2018b. Any suggestion on how to resolve this issue?

syms s K
assume(K,'real')
assumeAlso(K>0)
Sol=solve(K*s^3*(s+1)+1==0,s);

wrong_data=vpa(subs(Sol(1),K,linspace(10,20,100)));
right_data=double(subs(Sol(1),K,linspace(10,20,100)));

figure
subplot(2,1,1)
plot(real(right_data),imag(right_data),'bo');
subplot(2,1,2)
plot(real(wrong_data),imag(wrong_data),'r*');

[Walter Roberson]

I see what you mean. I will write an observation about this a little later this evening.

[S H]

Thanks. It is a week I am trying all possible approaches to get correct results without any success. Both VPA and double are erroneous on very large matrix data.

Answer 1

[Walter Roberson]

When you solve() symbolically you get the four roots symbolically,

root(K*z^4 + K*z^3 + 1, z, 1)
root(K*z^4 + K*z^3 + 1, z, 2)
root(K*z^4 + K*z^3 + 1, z, 3)
root(K*z^4 + K*z^3 + 1, z, 4)

where root(f(z), z) stands for the set of values, z, such that f(z) == 0 -- the roots of the equation.

The number after the z in root() is a root selector. For any given fully-numeric root, the symbolic toolbox is consistent about what value the root selector refers to.

Unfortunately, it is not documented as to which root a particular index corresponds to. root(f(x)+0,x,1) could potentially be a completely different branch than root(f(x)+epsilon,x,1) as far as the documentation is concerned.

What is known in practice about the root selectors? Not a lot, other than that the branch associated with any one index value is known to potentially change as some parameter changes.

syms z p
fplot(imag(root(z^4+z^3+p,z,1)),[-10 10])

Over the years I have noticed that for cubics, that if there is at least one real root, then root #1 appears to be real valued. But I do not know if that is a rule or if I have just not happened to try an exception.

Because of this, it would not at all surprise me if root(z^4+z^3+p,z,1) had discontinuities with changes in p (p = 1/K in your original notation.)

I have never before seen the fact that you point out, that the value obtained by applying double() to a particular root selector can be different than the value associated with vpa() of the same root selector, but my tests confirm your results. What you are observing is indeed present.

I can weakly point out that with the branch associated with a root index being undocumented, we can't really say that double() of a root() will be the same as vpa() of the root.

What I can point out a little more robustly is that vpa() operates at higher default digits, so algorithms that are in any way iterative could potentially end up on a different side of an ill-conditioned expression. I seem to remember that a couple of months ago I encountered a proof that root finding in polynomials is inherently tied to precision -- that there is no precision that is "good enough", that for the same polynomial, the root order will switch places an indefinite number of times as the precision increases; unfortunately at the moment I am a bit fuzzy on the details (and possibly the fundamental facts as well.)

So... what can you do? In this particular case you can do this:

Sol=solve(K*s^3*(s+1)+1==0,s, 'MaxDegree', 4);

then Sol(1) selects a particular branch in closed form.

These are pretty long expressions, and subs() of your linspace into them can take a long long time because it is going to try to find the algebraic number in each case. In practice you would matlabFunction() the formulae ... keeping in mind that there are a number of divisions in the formulae and that as you changed parameters you could encounter one of the singularities in lower precision that you would not necessarily encounter with higher precision.

From here we can ask what relationship we observe between the root order returned by double() of the root() forms, versuse vpa() or the root() forms, versus the direct closed form formulae. In the test values I ran through, the relationship was:

• vpa() order was exact reverse of the closed form order;
• double() order was iregular compared to closed form.

This relationship might not hold over other values or other polynomials.

For now I would say that if consistency is important to your purposes, then use MaxDegree to create the closed forms. This will not work for most polynomials of degree 5 or higher.

Over the longer term, I would suggest that you should consider selecting the root according to its properties. For example you might need the one with positive (non-zero) imaginary coefficient, or you might need the positive real-valued root. vpa() them all and filter afterwards.. and watch out for 0 matches and multiple matches.

[S H]

Thank you Walter for explaining this. To me it is still a huge bug that needs to be fixed for the following reasons.

The MaxDegree has an upper limit of 5. For polynomial orders higher than 5, the issue persists.

Matlab is doing something that is not clearly documented or explained to its customers. That can be traced by a simple root locus example. The rlocus performs similar function to solving equations as I have pointed out in my question above. The Matlab rlocus shows all results properly wothout any discontinuity in phase or imaginary part. But on the other hand, customers like me should solve our equations with the limitations that at different epsilons, the results might come from other root(z) solutions without any control.

So what is it that Matlab does in rlocus that avoids selecting datasets from different root(z) selections?

To make it more clear, see the plots from following code:

syms s K
tfunc=(s^2+2*s+2)/(s^5+9*s^4+33*s^3+51*s^2+26*s);
[symNum,symDen]=numden(tfunc);
TFnum=sym2poly(symNum);
TFden=sym2poly(symDen);
sys=tf(TFnum,TFden);
figure
rlocusplot(sys)
Sol=solve(K*symNum+symDen==0,s);
data=double(subs(Sol,K,logspace(-3,3,100)));
figure
hold on
for k=1:size(data,1)
    plot(real(data(k,:)),imag(data(k,:)))
end

Plotting root locus is very basic and its theory is known by general pubplic. What is wrong in above code that does not generate the correct plot? In another word, what is the undocumented trick that Matlab does in generating the correct results? Why results in Matlab rlocus plot does not come from different root(z) solutions.

[Walter Roberson]

I do not know much about the theory of control systems.

Tracing through to find the place the graph is computed was tricky. It turns out that it is not explicitly computed in an obvious way. Instead, rlocusplot sets a DataFcn property and when the graph is made visible a callback is made to the DataFcn which is then responsible for doing the real calculation.

The DataFcn that is set is toolbox/shared/controllib/graphics/@resppack/@ltisource/rlocus.m but it calls toolbox/control/ctrlobsolete/rlocus.m to do the real work -- Yes, even though it is marked as being in the "obsolete" section, it is used.

The "obsolete" rlocus routine calls ss() to convert num/denom form to state space form. That involves calling toolbox/shared/controllib/engine/numerics/compreal.m to calculate a number of companion matrices and balance them. I do not follow the mathematics of that part.

Then rlocus calls toolbox/shared/controllib/engine/+ltipack/@ssdata/utGetLoopData.m on the state space version. Part of the work involves,

% Reduce state-space data to Hessenberg form (stabilizes trajectories of 
% roots going to Inf)

Then the poles are calculated as the eigenvalues of one of the matrices, and ltipack.sszero calculates the zeros.

From there toolbox/shared/controllib/engine/+ltipack/@ltidata/private/gainrloc.m is responsible for computing the smooth graphs. It does some things I do not follow with gains and asymptopes . That calls the internal function smoothloc() which looks like it iterates through doing some kind of least squared matching to try to make the transitions smooth.

Thus, there is no internal use of solve() and double() or vpa() or any kind of closed form root formation that somehow results in is smooth following of root branches for rlocusplot(). The mathematics does appear to take some advantage of knowing what it computing in order to figure out which boundaries are important, but otherwise the magic comes down to generating roots and at boundaries selecting the smooth continuations adaptively. Figuring out the boundaries and understanding how the ss matrices work is significant for the methods, but the methods do not inherently follow branches: the code has to figure out how to follow the branches.

[S H]

Thank you walter for tracing the rlocus algorithm very clearly in detsails. I will look into those files you pointed out.

I also looked at selecting the roots based on their properties. For high order systems when there are many real or complex roots, separating them based on the sign of their values becomes impossible.

Again, thank you very much for taking the time and looking at the issue I raised.