How to reorganize a mathematical expression into its polynomial form

Viewed 67

I'm trying to reorganize the expression polin:

polin=1+K*(1+taua*s)/((1+tau1*s)*(1+tau2*s))*(1-s*theta/2)/(1+s*theta/2)*Kc*(1+1/(s*tauI)+tauD*s)

to this form: (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e)

I need to find the values of a, b, c, d, and e and it's complicated to do it by hand.

This is the code I used

import numpy as np
from sympy import *
Kc, K, tauI, taua, theta, tauD, tau1,tau2, s = sympy.symbols('Kc K tauI taua theta tauD tau1 tau2 s')

polin=1+K*(1+taua*s)/((1+tau1*s)*(1+tau2*s))*(1-s*theta/2)/(1+s*theta/2)*Kc*(1+1/(s*tauI)+tauD*s)

pol=poly_from_expr(polin,gen=s)
  
print(pol)

This is the output

(Poly(-1/2*s**3*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*taua*theta + s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*taua - 1/2*s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD*theta - 1/2*s**2*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*taua*theta - s*(1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI))*K*Kc*taua*theta + s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*tauD + s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*taua - 1/2*s*(1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc*theta - (1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI))*K*Kc*theta + (1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI))*K*Kc*taua + (1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1))*K*Kc + (1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI))*K*Kc + 1, s, 1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI), 1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI), 1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1), 1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI), K, Kc, tauD, taua, theta, domain='QQ'), {'gen': s, 'gens': (s, 1/(s**3*tau1*tau2*tauI*theta + 2*s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta + s**2*tau2*tauI*theta + 2*s*tau1*tauI + 2*s*tau2*tauI + s*tauI*theta + 2*tauI), 1/(s**3*tau1*tau2*tauI*theta/2 + s**2*tau1*tau2*tauI + s**2*tau1*tauI*theta/2 + s**2*tau2*tauI*theta/2 + s*tau1*tauI + s*tau2*tauI + s*tauI*theta/2 + tauI), 1/(s**3*tau1*tau2*theta/2 + s**2*tau1*tau2 + s**2*tau1*theta/2 + s**2*tau2*theta/2 + s*tau1 + s*tau2 + s*theta/2 + 1), 1/(s**4*tau1*tau2*tauI*theta/2 + s**3*tau1*tau2*tauI + s**3*tau1*tauI*theta/2 + s**3*tau2*tauI*theta/2 + s**2*tau1*tauI + s**2*tau2*tauI + s**2*tauI*theta/2 + s*tauI), K, Kc, tauD, taua, theta), 'domain': QQ, 'polys': False})

Thank you

1 Answers

Your expression looks like the following:

enter image description here

The order of terms in the polynomial in the numerator goes from -1 to 3 while the denominator goes from 0 to 3. This means that it is bottom-heavy and you will get a 1/s term. This means that there does not exist constants a, b, c, d, e such that polin == (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e).

The best you can do is an approximation using a series expansion. It is important to note that one of the terms is 1/s and so you will get something of the form (a)*s^4+ (b)*s^3+ (c)*s^2 + (d)*s + (e) + (f)/s.

To solve this, you could manually get the coefficients or multiply polin by s so that all powers are nonnegative.

Method 1:

from sympy import *

Kc, K, tauI, taua, theta, tauD, tau1, tau2, s = symbols('K_c K tau_I tau_a theta tau_D tau_1 tau_2 s', real=True)

polin = 1 + K * (1 + taua * s) / ((1 + tau1 * s) * (1 + tau2 * s)) * (1 - s * theta / 2) / (1 + s * theta / 2) * Kc * (
            1 + 1 / (s * tauI) + tauD * s)

n=5  # highest power of polin is s**(n-1)
x = series(polin, s, n=n).removeO()

coefficients = [x.coeff(s**(i)) for i in range(-1, n)]
for i, coef in enumerate(coefficients):
    print(f"Coefficient of s**{i-1} is {coef}")

Method 2:

from sympy import *

Kc, K, tauI, taua, theta, tauD, tau1, tau2, s = symbols('K_c K tau_I tau_a theta tau_D tau_1 tau_2 s', real=True)

polin = 1 + K * (1 + taua * s) / ((1 + tau1 * s) * (1 + tau2 * s)) * (1 - s * theta / 2) / (1 + s * theta / 2) * Kc * (
            1 + 1 / (s * tauI) + tauD * s)

n=6  # highest power of polin is s**(n-2)
x = series(s*polin, s, n=n).removeO()

sp: Poly = Poly(x, s)  # s multiplied by polin
coefficients = sp.all_coeffs()
for i, coef in enumerate(coefficients):
    print(f"Coefficient of s**{n-i-2} is {coef}")
Related