diff function rearranges symbolic expression leading to non-vanishing terms

I am trying to feedback linearize my non-linear system in Matlab. To that end I have to represent it in the following form:

So after I have modeled my system in the form of nonlinear differential equations, , I use the diff function to take the derrivative with respect to u to aquire g, and . This works, of course only if the system can be represented in the mentioned form. However, the diff function rewrites some of my expressions, which then after subtraction lead to non-vanishing factors multiplying u. For example, my input is still present in multiplied with a factor 5.7024e-18. I suspect, there is some quantization error introduced due to diff re-writing the derivative . Can this be somehow avoided?

g_x=[diff(dxdt,[Q_g]) diff(dxdt,[Q_b])];
for i=1:size(g_x,1)
    i
f_x(i,:)=simplify(dxdt(i,:)-g_x(i,1)*Q_g-g_x(i,2)*Q_b,'steps',10,'IgnoreAnalyticConstraints',true);

end   

Here Q_g and Q_b are my inputs.

An example of the re-writing of expressions:

dxdt = - (1297036692682702848*Q_b*(C_carbv - 47/20000))/15694143981460705 + ..................

g_x=1905022642377719808/9808839988412940625 - (1297036692682702848*C_carbv)/15694143981460705

f_x=Q_b*((1297036692682702848*C_carbv)/15694143981460705 - 1905022642377719808/9808839988412940625) - (1297036692682702848*Q_b*(C_carbv - 47/20000))/15694143981460705 + ..............................

simplify(Q_b*((1297036692682702848*C_carbv)/15694143981460705 - 1905022642377719808/9808839988412940625) - ...
    (1297036692682702848*Q_b*(C_carbv - 47/20000))/15694143981460705)

ans =

(47*Q_b)/8242150268041428750

The .................. stands for terms that do not depend on Q_b.

[Paul]

Can you provide complete code that can be replicated that illustrates the problem and is prefably as simplie as possible?

[Martin Elenkov]

This is a script with the model parameters and equations. The cited expression are under index 9. So, dxdt(9), f_x(9), g_x(9,1:2). It is not simple, nor clean, but I think it is fairly clear what I am doing, asside from the complicated differential equations

Answer 1

[Paul]

Why not just use collect() and coeffs() on the elements of dxdt? If the system is affine in u, then those should give the desired result directly

syms x u1 u2
dxdt = x^2 + x + x*u1 + x^3*u1 + x^2*u2 + u2;
dxdt = collect(dxdt,[u1 u2])

dxdt = 

[c,t] = coeffs(dxdt,[u1 u2])

c = 

t = 

fofx = c(find(t==1))

fofx = 

gofx = [c(find(t==u1)) c(find(t==u2))]

gofx = 

If the system is not affine in u, then extra terms will show up in t, so you can check numel(t) before doing any other processing to decide how to proceed.

dxdt = x^2 + x + x*u1 + x^3*u1 + x^2*u2 + u2 + u1*u2;
dxdt = collect(dxdt,[u1 u2])

dxdt = 

[c,t] = coeffs(dxdt,[u1 u2])

c = 

t = 

[Martin Elenkov]

So elegant! So easy! Thank you Paul!

[Paul]

collect() was just for illustration. Only coeffs() is needed.

[Walter Roberson]

If you use the 'all' option of coeffs, then t = 1 will always be the last coefficient

[Paul]

I thought 'all' might also make it easier to determine if f(x) or an element of g(x) is zero. Alas, that appears to not be the case:

syms x u1 u2
dxdt = x*u1;
[c,t]=coeffs(dxdt,[u1 u2],'all')

c = 

t = 

So some logic will still be needed in a situation like this to set the appropriate element of g(x) to zero because u2 doesn't exist in t, even though it was explicitly requested, which is too bad for this use case.

[Martin Elenkov]

I also noticed that only the existing elements are returned, but it is an easy work around.

Answer 2

[John D'Errico]

This is a floating point problem, but it comes when MATLAB takes your constants and turns them into symbolic numbers. Do you see all of those large integers? For example...

X = sym(1.34536363425474)

X = 

But that ratio is not exactly the same value as the original number. It is typically a close approximation.

vpa(X,40)

ans = 

1.345363634254739926277011363708879798651

So you get trash in there, down in digit 17 and below.

My question is, why is this a problem?

[Walter Roberson]

When you have a floating point value involved in a symbolic expression, then by default MATLAB converts it using 'r' option.

The 'r' option asks the symbolic toolbox to take the floating point number and run a series of approximations and return one of:

• a simple rational multiple of πwith a denominator no greater than 100000 (most larger denominators will not be detected)
• a simple rational multiple of a sqrt() of a rational integer; when the denominator starts getting larger than about 315 then it starts occasionally missing the possibility that the number can be expressed this way
• or, if need be, a rational number

In each case, numbers that are very close to the above will be approximated as being exactly one of the above. For example sym(pi*(1-100*eps)) will be approximated as πexactly but sym(pi*(1-101*eps)) will not be

If you have a general number that has no obvious structure, such as the 1.34536363425474 that John mentions, then most of the time it is going to be expressed as a rational with large denominator.

You can, however, use sym() yourself to ask to transform a floating point number to symbolic in a different way. The most common way to do that is to use

sym(NUMBER, 'd')

This does not do any approximation. Instead, it takes the exact floating point binary representation, and converts it into decimal according to the current digits setting -- so for example sym(1/6,'d') is not going to show up as 1.666 with the 6 repeated until the end, but will instead show up like 0.1666666666666666574148081281237 because 1/6 had to be approximated to store as floating point, and it is the approximated value that is converted.

You can also use sym(NUMBER, 'f')

This does not do any approximation: instead it takes the floating point representation in binary and expresses it as a rational, like 6004799503160661/36028797018963968 ... remembering that 1/6 has no exact finite binary representation

Sometimes what you end up wanting to do is

vpa(sym(NUMBER))

so for example sym(1/6) invokes the default 'r' rational approximation to find the symbolic rational number 1/6 and then vpa() can approximate that into decimal for as many decimal places as needed. Likewise vpa(sym(1.57)) has the sym(1.57) find the rational approximation 157/100 and then vpa() of that expresses it as symbolic floating point 1.57

And sometimes it is best to just sym() of a quoted number, such as

sym('1.34536363425474')

to get the same decimal number out as software floating point.

[Martin Elenkov]

@John D'Errico The problem is that if I solve this analytically I will aquire as a function of x only, but matlab returns a function , which also depends on the input. I started with a function of the states and inputs and want to rewrite it into one function depending only on the states and one that depends on the states and the inputs. It seems like I have to treat the constants in a different manner - for example, approximate them with vpa()

@Walter Roberson Thanks for your thorough answer! I tried:

f_x=simplify(vpa(dxdt,4)-vpa(g_x,4)*Q_g-vpa(g_x,4)*Q_b)

which unexpectedly returns numbers that have much more precision than 4 digits:

Q_b*(82.644628099165856838226318359375*C_carbv - 0.19421487603312925784848630428314) - 82.644628099165856838226318359375*Q_b*(C_carbv - 0.0023499999999998522071109618991613)

vpa(f_x,2)

Q_b*(83.0*C_carbv - 0.19) - 83.0*Q_b*(C_carbv - 2.3e-3)

In the last result, you can clearly see that the expressiosns should be equal, however MATLAB still does not want to do the cancelation opperation.

Should I be using vpa() much earlier? i.e defining the model parameters as symbolic variables?

[Walter Roberson]

simplify() asks to rewrite the expression to combine parts if it can figure out how. When it does that, the combined portions mostly lose the information about the number of digits the symbolic floating point numbers were evaluated at

syms Q_b C_carbv 
f_x = Q_b * (vpa(83.0,2)*C_carbv - vpa(0.19,2)) - vpa(83.0,2)*Q_b*(C_carbv - vpa(2.3e-3,2))

f_x = 

expand(f_x)

ans = 

Makes sense to me, since - 0.19 is not - 2.3e-3

collect(f_x, Q_b)

ans = 

[Martin Elenkov]

Notice that in the first term Matlab factors only in front of the brackets, while in the second term, both the coefficient 83 and are factored out of the brackets (They both come from the same equation, Matlab just rewrites and rearranges them differently after some opperations - e.g. diff()).

83*2.3e-3

ans =

0.1909

0.19 is not 2.3e-3, but 83*2.3e-3 is aproximately equal to 0.19, depending on how you define your precision. It seems that Matlab reaproximates these coefficient after each (or at least some) opperations and that leads to this anoying problem.