Evaluate how we update the speed during numerical integration

[raphaeltimbo]

Our space state matrix A has to be assembled every time we change the speed. This can be a problem for numerical integration.

[ross-bott]

Hi there! I have marked this issue as stale because it has not had activity for 45 days. Consider the following options:

• If the issue refers to a large task, break it in smaller issues that can be solved in less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.

[Raimundovpn]

The matrix of state spaces “A” depends on the speed and acceleration of rotation, however, in the current method, for each speed of rotation the matrix “A” is reassembled and reintegrated. To avoid this, it would be necessary to transform the result of the integration (the time response) into a function of the speed and acceleration of rotation. The state space matrix should be reformulated, in order to isolate the terms that depend on the speed and acceleration of the rotor, and consequently, modify the integration method, as it currently uses scipy. The result of this would be a function that could be evaluated at any speed/acceleration of rotation without the need to reintegrate matrix A.

[Raimundovpn]

The state space matrix defined by the function def A(self, speed=0, frequency=None): depends on the speed and acceleration of rotation, however, by the current method, it needs to be assembled and integrated for each speed of rotation , because the speed is a parameter of the A matrix.

def A(self, speed=0, frequency=None):
        """State space matrix for an instance of a rotor.

        Parameters
        ----------
        speed: float, optional
            Rotor speed.
            Default is 0.
        frequency : float, optional
            Excitation frequency. Default is rotor speed.

        Returns
        -------
        A : np.ndarray
            State space matrix for the rotor.

        Examples
        --------
        >>> rotor = rotor_example()
        >>> np.round(rotor.A()[50:56, :2])
        array([[     0.,  10927.],
               [-10924.,     -0.],
               [  -174.,      0.],
               [    -0.,   -174.],
               [    -0.,  10723.],
               [-10719.,     -0.]])
        """
        if frequency is None:
            frequency = speed

        Z = np.zeros((self.ndof, self.ndof))
        I = np.eye(self.ndof)

        # fmt: off
        A = np.vstack(
            [np.hstack([Z, I]),
             np.hstack([la.solve(-self.M(), self.K(frequency) + self.Kst()*speed), la.solve(-self.M(), (self.C(frequency) + self.G() * speed))])])
        # fmt: on

        return A

According to the current method, the integration process first consists of building a linear system in time, using the scipy.signal library and from there, determining the other analyses. The method is defined as follows:

 def _lti(self, speed, frequency=None):
        """Continuous-time linear time invariant system.

        This method is used to create a Continuous-time linear
        time invariant system for the mdof system.
        From this system we can obtain poles, impulse response,
        generate a bode, etc.

        Parameters
        ----------
        speed: float
            Rotor speed.
        frequency: float, optional
            Excitation frequency.
            Default is rotor speed.

        Returns
        -------
        sys : StateSpaceContinuous
            Space State Continuos with A, B, C and D matrices

        Example
        -------
        >>> rotor = rotor_example()
        >>> A = rotor._lti(speed=0).A
        >>> B = rotor._lti(speed=0).B
        >>> C = rotor._lti(speed=0).C
        >>> D = rotor._lti(speed=0).D
        """
        Z = np.zeros((self.ndof, self.ndof))
        I = np.eye(self.ndof)

        # x' = Ax + Bu
        B2 = I
        if frequency is None:
            frequency = speed
        A = self.A(speed=speed, frequency=frequency)
        # fmt: off
        B = np.vstack([Z,
                       la.solve(self.M(), B2)])
        # fmt: on

        # y = Cx + Du
        # Observation matrices
        Cd = I
        Cv = Z
        Ca = Z

        # fmt: off
        C = np.hstack((Cd - Ca @ la.solve(self.M(), self.K(frequency)), Cv - Ca @ la.solve(self.M(), self.C(frequency))))
        # fmt: on
        D = Ca @ la.solve(self.M(), B2)

        sys = signal.lti(A, B, C, D)

        return sys

In addition, this function is used to build the functions: def time_response(), def transfer_matrix(), and through them to carry out the other analyses.

The main advantage of the current method is its simplicity. This makes the code very clear and easy to understand, the problem is that it is not very flexible, which makes it impossible to take advantage of some structures already calculated in other analyses. For example, in the state space matrix def A(), the matrices that compose it do not change with a new rotation speed, the ideal would be to rewrite it, in order to isolate the terms that depend on the speed and acceleration of the rotor, and consequently modify the integration method, as it currently uses scipy. The result of this would be a function that could be evaluated at any speed/acceleration of rotation without the need to reintegrate matrix A, thus reducing the time and computational cost of each analysis that depended on it.

[raphaeltimbo]

Closing this since we already have #1021 to discuss this topic.