AndreyBychkov / QBee

Quadratization of differential equations
MIT License
7 stars 2 forks source link

Received message -> RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode. #23

Open smallpondtom opened 4 months ago

smallpondtom commented 4 months ago

Hi. I'm trying out a pretty ambitious ODE from orbital mechanics called the Modified Equinoctial Element experimentally. But I received the message

RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode.

So, I'm letting you guys know.

You can reproduce this with the following code:

import numpy as np
import sympy as sp
from qbee import * 

# Define the parameters
mu = parameters("mu")

# Define the states 
p, f, g, h, k, L, D1, D2, D3 = functions("p f g h k L D1 D2 D3")

# Define additional auxiliary variables
q = 1 + f*sp.cos(L) + g*sp.sin(L)
s2 = 1 + h**2 + k**2

# Define the Modified Equinoctial Elements
pdot = 2*p*sp.sqrt(p/mu)/q * D2
fdot = sp.sqrt(p/mu)*sp.sin(L)*D1 + sp.sqrt(p/mu)/q*((q+1)*sp.cos(L) + f)*D2 - sp.sqrt(p/mu)*g/q*(h*sp.cos(L) - k*sp.sin(L))*D3
gdot = -sp.sqrt(p/mu)*sp.cos(L)*D1 + sp.sqrt(p/mu)/q*((q+1)*sp.sin(L) + g)*D2 + sp.sqrt(p/mu)*f/q*(h*sp.cos(L) - k*sp.sin(L))*D3
hdot = sp.sqrt(p/mu)*s2*sp.cos(L)/2/q * D3
kdot = sp.sqrt(p/mu)*s2*sp.sin(L)/2/q * D3
Ldot = sp.sqrt(p/mu)/q*(h*sp.sin(L) - k*sp.cos(L)) * D3
b6 = sp.sqrt(mu*p)*q**2/p**2 

system = [
    (p, pdot),
    (f, fdot),
    (g, gdot),
    (h, hdot),
    (k, kdot),
    (L, Ldot + b6)
]

polynomialize_and_quadratize(system, input_free=True).print()

Thanks.

pogudingleb commented 4 months ago

@smallpondtom Thanks for reporting this ! Just to clarify that I understand the system correctly: D1, D2, D3 are external time-varying inputs, is that correct?

smallpondtom commented 4 months ago

@pogudingleb Yes. I'm not an expert on orbital mechanics but I think the control inputs are essentially like thruster/actuator inputs. And, they are time-varying but not necessarily continuous since it could be a simple on/off control. That's why I used the input_free option.

pogudingleb commented 4 months ago

@smallpondtom Thanks again for the example! There is a bug in the polynomialization part which we will fix. However, even if I help the software to get polynomialization right, the model is too large and too nonlinear for our quadratization algorithm.