AutoAndroidController/python/Lib/site-packages/sympy/physics/control/tests/test_control_plots.py

from math import isclose
from sympy.core.numbers import I, all_close
from sympy.core.symbol import Dummy
from sympy.functions.elementary.complexes import (Abs, arg)
from sympy.functions.elementary.exponential import log
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.abc import s, p, a
from sympy import pi
from sympy.external import import_module
from sympy.physics.control.control_plots import \
    (pole_zero_numerical_data, pole_zero_plot, step_response_numerical_data,
    step_response_plot, impulse_response_numerical_data,
    impulse_response_plot, ramp_response_numerical_data,
    ramp_response_plot, bode_magnitude_numerical_data,
    bode_phase_numerical_data, bode_plot, nyquist_plot_expr,
    nichols_plot_expr)

from sympy.physics.control.lti import (TransferFunction,
    Series, Parallel, TransferFunctionMatrix)
from sympy.testing.pytest import raises, skip

matplotlib = import_module(
        'matplotlib', import_kwargs={'fromlist': ['pyplot']},
        catch=(RuntimeError,))

numpy = import_module('numpy')

tf1 = TransferFunction(1, p**2 + 0.5*p + 2, p)
tf2 = TransferFunction(p, 6*p**2 + 3*p + 1, p)
tf3 = TransferFunction(p, p**3 - 1, p)
tf4 = TransferFunction(10, p**3, p)
tf5 = TransferFunction(5, s**2 + 2*s + 10, s)
tf6 = TransferFunction(1, 1, s)
tf7 = TransferFunction(4*s*3 + 9*s**2 + 0.1*s + 11, 8*s**6 + 9*s**4 + 11, s)
tf8 = TransferFunction(5, s**2 + (2+I)*s + 10, s)

ser1 = Series(tf4, TransferFunction(1, p - 5, p))
ser2 = Series(tf3, TransferFunction(p, p + 2, p))

par1 = Parallel(tf1, tf2)


def _to_tuple(a, b):
    return tuple(a), tuple(b)

def _trim_tuple(a, b):
    a, b = _to_tuple(a, b)
    return tuple(a[0: 2] + a[len(a)//2 : len(a)//2 + 1] + a[-2:]), \
        tuple(b[0: 2] + b[len(b)//2 : len(b)//2 + 1] + b[-2:])

def y_coordinate_equality(plot_data_func, evalf_func, system):
    """Checks whether the y-coordinate value of the plotted
    data point is equal to the value of the function at a
    particular x."""
    x, y = plot_data_func(system)
    x, y = _trim_tuple(x, y)
    y_exp = tuple(evalf_func(system, x_i) for x_i in x)
    return all(Abs(y_exp_i - y_i) < 1e-8 for y_exp_i, y_i in zip(y_exp, y))


def test_errors():
    if not matplotlib:
        skip("Matplotlib not the default backend")

    # Invalid `system` check
    tfm = TransferFunctionMatrix([[tf6, tf5], [tf5, tf6]])
    expr = 1/(s**2 - 1)
    raises(NotImplementedError, lambda: pole_zero_plot(tfm))
    raises(NotImplementedError, lambda: pole_zero_numerical_data(expr))
    raises(NotImplementedError, lambda: impulse_response_plot(expr))
    raises(NotImplementedError, lambda: impulse_response_numerical_data(tfm))
    raises(NotImplementedError, lambda: step_response_plot(tfm))
    raises(NotImplementedError, lambda: step_response_numerical_data(expr))
    raises(NotImplementedError, lambda: ramp_response_plot(expr))
    raises(NotImplementedError, lambda: ramp_response_numerical_data(tfm))
    raises(NotImplementedError, lambda: bode_plot(tfm))

    # More than 1 variables
    tf_a = TransferFunction(a, s + 1, s)
    raises(ValueError, lambda: pole_zero_plot(tf_a))
    raises(ValueError, lambda: pole_zero_numerical_data(tf_a))
    raises(ValueError, lambda: impulse_response_plot(tf_a))
    raises(ValueError, lambda: impulse_response_numerical_data(tf_a))
    raises(ValueError, lambda: step_response_plot(tf_a))
    raises(ValueError, lambda: step_response_numerical_data(tf_a))
    raises(ValueError, lambda: ramp_response_plot(tf_a))
    raises(ValueError, lambda: ramp_response_numerical_data(tf_a))
    raises(ValueError, lambda: bode_plot(tf_a))

    # lower_limit > 0 for response plots
    raises(ValueError, lambda: impulse_response_plot(tf1, lower_limit=-1))
    raises(ValueError, lambda: step_response_plot(tf1, lower_limit=-0.1))
    raises(ValueError, lambda: ramp_response_plot(tf1, lower_limit=-4/3))

    # slope in ramp_response_plot() is negative
    raises(ValueError, lambda: ramp_response_plot(tf1, slope=-0.1))

    # incorrect frequency or phase unit
    raises(ValueError, lambda: bode_plot(tf1,freq_unit = 'hz'))
    raises(ValueError, lambda: bode_plot(tf1,phase_unit = 'degree'))


def test_pole_zero():

    def pz_tester(sys, expected_value):
        _z, _p = pole_zero_numerical_data(sys)
        z_check = all_close(_z, expected_value[0])
        p_check = all_close(_p, expected_value[1])
        return p_check and z_check

    exp1 = [[], [-0.24999999999999994-1.3919410907075054j, -0.24999999999999994+1.3919410907075054j]]
    exp2 = [[0.0], [-0.25-0.3227486121839514j, -0.25+0.3227486121839514j]]
    exp3 = [[0.0], [0.9999999999999998+0j, -0.5000000000000004-0.8660254037844395j,
        -0.5000000000000004+0.8660254037844395j]]
    exp4 = [[], [0.0, 0.0, 0.0, 5.0]]
    exp5 = [[-5.645751311064592, -0.5000000000000008, -0.3542486889354093],
        [-0.24999999999999986-0.322748612183951348j,
        -0.2499999999999998+0.32274861218395134j,
        -0.24999999999999986-1.3919410907075052j,
         -0.2499999999999998+1.3919410907075052j]]
    exp6 = [[], [-1.1641600331447917-3.545808351896439j,
          -0.8358399668552097+2.5458083518964383j]]

    assert pz_tester(tf1, exp1)
    assert pz_tester(tf2, exp2)
    assert pz_tester(tf3, exp3)
    assert pz_tester(ser1, exp4)
    assert pz_tester(par1, exp5)
    assert pz_tester(tf8, exp6)


def test_bode():
    if not numpy:
        skip("NumPy is required for this test")

    def bode_phase_evalf(system, point):
        expr = system.to_expr()
        _w = Dummy("w", real=True)
        w_expr = expr.subs({system.var: I*_w})
        return arg(w_expr).subs({_w: point}).evalf()

    def bode_mag_evalf(system, point):
        expr = system.to_expr()
        _w = Dummy("w", real=True)
        w_expr = expr.subs({system.var: I*_w})
        return 20*log(Abs(w_expr), 10).subs({_w: point}).evalf()

    def test_bode_data(sys):
        return y_coordinate_equality(bode_magnitude_numerical_data, bode_mag_evalf, sys) \
            and y_coordinate_equality(bode_phase_numerical_data, bode_phase_evalf, sys)

    assert test_bode_data(tf1)
    assert test_bode_data(tf2)
    assert test_bode_data(tf3)
    assert test_bode_data(tf4)
    assert test_bode_data(tf5)


def check_point_accuracy(a, b):
    return all(isclose(*_, rel_tol=1e-1, abs_tol=1e-6
        ) for _ in zip(a, b))


def test_impulse_response():
    if not numpy:
        skip("NumPy is required for this test")

    def impulse_res_tester(sys, expected_value):
        x, y = _to_tuple(*impulse_response_numerical_data(sys,
            adaptive=False, n=10))
        x_check = check_point_accuracy(x, expected_value[0])
        y_check = check_point_accuracy(y, expected_value[1])
        return x_check and y_check

    exp1 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (0.0, 0.544019738507865, 0.01993849743234938, -0.31140243360893216, -0.022852779906491996, 0.1778306498155759,
        0.01962941084328499, -0.1013115194573652, -0.014975541213105696, 0.0575789724730714))
    exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.1666666675, 0.08389223412935855,
        0.02338051973475047, -0.014966807776379383, -0.034645954223054234, -0.040560075735512804,
        -0.037658628907103885, -0.030149507719590022, -0.021162090730736834, -0.012721292737437523))
    exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (4.369893391586999e-09, 1.1750333000630964,
        3.2922404058312473, 9.432290008148343, 28.37098083007151, 86.18577464367974, 261.90356653762115,
        795.6538758627842, 2416.9920942096983, 7342.159505206647))
    exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 6.17283950617284, 24.69135802469136,
        55.555555555555564, 98.76543209876544, 154.320987654321, 222.22222222222226, 302.46913580246917,
        395.0617283950618, 500.0))
    exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, -0.10455606138085417,
        0.06757671513476461, -0.03234567568833768, 0.013582514927757873, -0.005273419510705473,
        0.0019364083003354075, -0.000680070134067832, 0.00022969845960406913, -7.476094359583917e-05))
    exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (-6.016699583000218e-09, 0.35039802056107394, 3.3728423827689884, 12.119846079276684,
        25.86101014293389, 29.352480635282088, -30.49475907497664, -273.8717189554019, -863.2381702029659,
        -1747.0262164682233))
    exp7 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335,
        4.444444444444445, 5.555555555555555, 6.666666666666667, 7.777777777777779,
        8.88888888888889, 10.0), (0.0, 18.934638095560974, 5346.93244680907, 1384609.8718249386,
        358161126.65801865, 92645770015.70108, 23964739753087.42, 6198974342083139.0, 1.603492601616059e+18,
        4.147764422869658e+20))

    assert impulse_res_tester(tf1, exp1)
    assert impulse_res_tester(tf2, exp2)
    assert impulse_res_tester(tf3, exp3)
    assert impulse_res_tester(tf4, exp4)
    assert impulse_res_tester(tf5, exp5)
    assert impulse_res_tester(tf7, exp6)
    assert impulse_res_tester(ser1, exp7)


def test_step_response():
    if not numpy:
        skip("NumPy is required for this test")

    def step_res_tester(sys, expected_value):
        x, y = _to_tuple(*step_response_numerical_data(sys,
            adaptive=False, n=10))
        x_check = check_point_accuracy(x, expected_value[0])
        y_check = check_point_accuracy(y, expected_value[1])
        return x_check and y_check

    exp1 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (-1.9193285738516863e-08, 0.42283495488246126, 0.7840485977945262, 0.5546841805655717,
        0.33903033806932087, 0.4627251747410237, 0.5909907598988051, 0.5247213989553071,
        0.4486997874319281, 0.4839358435839171))
    exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (0.0, 0.13728409095645816, 0.19474559355325086, 0.1974909129243011, 0.16841657696573073,
        0.12559777736159378, 0.08153828016664713, 0.04360471317348958, 0.015072994568868221,
        -0.003636420058445484))
    exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (0.0, 0.6314542141914303, 2.9356520038101035, 9.37731009663807, 28.452300356688376,
        86.25721933273988, 261.9236645044672, 795.6435410577224, 2416.9786984578764, 7342.154119725917))
    exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (0.0, 2.286236899862826, 18.28989519890261, 61.72839629629631, 146.31916159122088, 285.7796124828532,
        493.8271703703705, 784.1792566529494, 1170.553292729767, 1666.6667))
    exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (-3.999999997894577e-09, 0.6720357068882895, 0.4429938256137113, 0.5182010838004518,
        0.4944139147159695, 0.5016379853883338, 0.4995466896527733, 0.5001154784851325,
        0.49997448824584123, 0.5000039745919259))
    exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (-1.5433688493882158e-09, 0.3428705539937336, 1.1253619102202777, 3.1849962651016517,
        9.47532757182671, 28.727231099148135, 87.29426924860557, 265.2138681048606, 805.6636260007757,
        2447.387582370878))

    assert step_res_tester(tf1, exp1)
    assert step_res_tester(tf2, exp2)
    assert step_res_tester(tf3, exp3)
    assert step_res_tester(tf4, exp4)
    assert step_res_tester(tf5, exp5)
    assert step_res_tester(ser2, exp6)


def test_ramp_response():
    if not numpy:
        skip("NumPy is required for this test")

    def ramp_res_tester(sys, num_points, expected_value, slope=1):
        x, y = _to_tuple(*ramp_response_numerical_data(sys,
            slope=slope, adaptive=False, n=num_points))
        x_check = check_point_accuracy(x, expected_value[0])
        y_check = check_point_accuracy(y, expected_value[1])
        return x_check and y_check

    exp1 = ((0.0, 2.0, 4.0, 6.0, 8.0, 10.0), (0.0, 0.7324667795033895, 1.9909720978650398,
        2.7956587704217783, 3.9224897567931514, 4.85022655284895))
    exp2 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445,
        5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0),
        (2.4360213402019326e-08, 0.10175320182493253, 0.33057612497658406, 0.5967937263298935,
        0.8431511866718248, 1.0398805391471613, 1.1776043125035738, 1.2600994825747305, 1.2981042689274653,
        1.304684417610106))
    exp3 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (-3.9329040468771836e-08,
        0.34686634635794555, 2.9998828170537903, 12.33303690737476, 40.993913948137795, 127.84145222317912,
        391.41713691996, 1192.0006858708389, 3623.9808672503405, 11011.728034546572))
    exp4 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 1.9051973784484078, 30.483158055174524,
        154.32098765432104, 487.7305288827924, 1190.7483615302544, 2469.1358024691367, 4574.3789056546275,
        7803.688462124678, 12500.0))
    exp5 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 3.8844361856975635, 9.141792069209865,
        14.096349157657231, 19.09783068994694, 24.10179770390321, 29.09907319114121, 34.10040420185154,
        39.09983919254265, 44.10006013058409))
    exp6 = ((0.0, 1.1111111111111112, 2.2222222222222223, 3.3333333333333335, 4.444444444444445, 5.555555555555555,
        6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0), (0.0, 1.1111111111111112, 2.2222222222222223,
        3.3333333333333335, 4.444444444444445, 5.555555555555555, 6.666666666666667, 7.777777777777779, 8.88888888888889, 10.0))

    assert ramp_res_tester(tf1, 6, exp1)
    assert ramp_res_tester(tf2, 10, exp2, 1.2)
    assert ramp_res_tester(tf3, 10, exp3, 1.5)
    assert ramp_res_tester(tf4, 10, exp4, 3)
    assert ramp_res_tester(tf5, 10, exp5, 9)
    assert ramp_res_tester(tf6, 10, exp6)


def test_nyquist_plot_expr():
    r1, i1, w1 = nyquist_plot_expr(tf1)
    r2, i2, w2 = nyquist_plot_expr(tf2)
    r3, i3, w3 = nyquist_plot_expr(tf3)
    r4, i4, w4 = nyquist_plot_expr(tf4)
    assert r1 == (2 - w1**2)/(0.25*w1**2 + (2 - w1**2)**2)
    assert i1 == -0.5*w1/(0.25*w1**2 + (2 - w1**2)**2)
    assert r2 == 3*w2**2/(9*w2**2 + (1 - 6*w2**2)**2)
    assert i2 == w2*(1 - 6*w2**2)/(9*w2**2 + (1 - 6*w2**2)**2)
    assert r3 == -w3**4/(w3**6 + 1)
    assert i3 == -w3/(w3**6 + 1)
    assert r4 == 0
    assert i4 == 10/w4**3


def test_nichols_expr():
    m1, p1, w1 = nichols_plot_expr(tf1)
    m2, p2, w2 = nichols_plot_expr(tf2)
    m3, p3, w3 = nichols_plot_expr(tf3)
    m4, p4, w4 = nichols_plot_expr(tf4)
    assert m1 == 20*log(1/sqrt(w1**4 - 3.75*w1**2 + 4))/log(10)
    assert p1 == 180*arg(1/(-w1**2 + 0.5*w1*I + 2))/pi
    assert m2 == 20*log(Abs(w2)/sqrt(36*w2**4 - 3*w2**2 + 1))/log(10)
    assert p2 == 180*arg(w2*I/(-6*w2**2 + 3*w2*I + 1))/pi
    assert m3 == 20*log(Abs(w3)/sqrt(w3**6 + 1))/log(10)
    assert p3 == 180*arg(-w3*I/(w3**3*I + 1))/pi
    assert m4 == 20*log(10/(w4**2*Abs(w4)))/log(10)
    assert p4 == 180*arg(I/w4**3)/pi