from math import isclose
from sympy import lambdify, latex
from sympy.core.numbers import I, pi
from sympy.core.symbol import Dummy
from sympy.functions.elementary.complexes import (Abs, arg)
from sympy.functions.elementary.exponential import log
from sympy.abc import s, p, a
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)
from sympy.physics.control.lti import (TransferFunction,
Series, Parallel, TransferFunctionMatrix)
from sympy.polys.polytools import Poly
from sympy.plotting.series import LineOver1DRangeSeries
matplotlib = import_module(
'matplotlib', import_kwargs={'fromlist': ['pyplot']},
catch=(RuntimeError,))
numpy = import_module('numpy')
if matplotlib:
plt = matplotlib.pyplot
if numpy:
np = numpy # Matplotlib already has numpy as a compulsory dependency. No need to install it separately.
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)
x, magnitude_data = bode_magnitude_numerical_data(tf1)
rad, phase_data = bode_phase_numerical_data(tf1, phase_unit='deg')
print(len(magnitude_data[0: len(phase_data)]))
print(len(phase_data))
magnitude_data = magnitude_data[0: len(phase_data)] #cut some of the last data of magnitude because len(phase_data) = len(magnitude_data)
#plt.plot(phase_data, magnitude_data)
#plt.yscale('log')
#plt.show()
def nichols_numerical_data(system, initial_omega=0.01, final_omega=100, **kwargs):
"""
Returns the numerical data of Nichols plot of the system.
It is internally used by ``nichols_plot`` to get the data
for plotting Nichols plot. Users can use this data to further
analyse the dynamics of the system or plot using a different
backend/plotting-module.
Parameters
==========
system : SISOLinearTimeInvariant
The system for which the pole-zero data is to be computed.
initial_omega : Number, optional
The initial value of frequency. Defaults to 0.01.
final_omega : Number, optional
The final value of frequency. Defaults to 100.
Returns
=======
tuple : (phase_points, mag_points)
phase_points = phase values of the Nichols plot.
mag_points = magnitude values of the Nichols plot.
Raises
======
NotImplementedError
When a SISO LTI system is not passed.
When time delay terms are present in the system.
ValueError
When more than one free symbol is present in the system.
The only variable in the transfer function should be
the variable of the Laplace transform.
Examples
========
>>> from sympy.abc import s
>>> from sympy.physics.control.lti import TransferFunction
>>> from sympy.physics.control.control_plots import nichols_numerical_data
>>> tf1 = TransferFunction(-(0.1)*s**3 - (2.4)*s**2 - 181*s - 1950, s**3 + (3.3)*s**2 + 990*s + 2600, s)
>>> nichols_numerical_data(tf1) # doctest: +SKIP
(array([179.83501857, 179.67004337, 179.50508061, ..., 166.86071969,
166.86233751, 166.8639549 ]),
array([ -2.49883392, -2.49901149, -2.49930742, ..., -20.5300856,
-20.52996573, -20.52984591]))
See Also
========
nichols_plot
"""
#_check_system(system)
expr = system.to_expr()
_w = Dummy("w", real=True)
repl = I*_w
w_expr = expr.subs({system.var: repl})
mag = 20*log(Abs(w_expr), 10)
phase = arg(w_expr)*180/pi
x = np.linspace(initial_omega, final_omega, 10000)
#mag_func = lambdify(_w, mag)
#phase_func = lambdify(_w, phase)
hz, mag_points = LineOver1DRangeSeries(mag, x).get_points()
rad, phase_points = LineOver1DRangeSeries(phase, x).get_points()
mag_points = mag_points[0: len(phase_points)] #cut some of the last data of magnitude because len(phase_data) = len(magnitude_data)
return phase_points, mag_points
def nichols_plot(system, initial_omega=0.01, final_omega=100,
color='b', grid=False, show=True,**kwargs):
r"""
Returns the nichols plot of a continuous-time system.
Nichols Plot is a plot used in signal processing and control system design
to determine the stability of a feedback system
Parameters
==========
system : SISOLinearTimeInvariant type
The LTI SISO system for which the Ramp Response is to be computed.
initial_omega : Number, optional
The initial value of frequency. Defaults to 0.01.
final_omega : Number, optional
The final value of frequency. Defaults to 100.
show : boolean, optional
If ``True``, the plot will be displayed otherwise
the equivalent matplotlib ``plot`` object will be returned.
Defaults to True.
show_axes : boolean, optional
If ``True``, the coordinate axes will be shown. Defaults to False.
grid : boolean, optional
If ``True``, the plot will have a grid. Defaults to False.
Examples
========
.. plot::
:context: close-figs
:format: doctest
:include-source: True
>>> from sympy.abc import s
>>> from sympy.physics.control.lti import TransferFunction
>>> from sympy.physics.control.control_plots import nichols_plot
>>> tf1 = TransferFunction(-(0.1)*s**3 - (2.4)*s**2 - 181*s - 1950, s**3 + (3.3)*s**2 + 990*s + 2600, s)
>>> nichols_plot(tf1, 1, 100) # doctest: +SKIP
See Also
========
bode_magnitude_plot, bode_phase_plot
References
==========
"""
x, y = nichols_numerical_data(system, initial_omega=initial_omega,
final_omega=final_omega)
plt.plot(x, y, color=color, **kwargs)
plt.xlabel('Open Loop Gain (dB)')
plt.ylabel('Open Loop Phase (deg)')
plt.title(f'Nichols Plot (Phase) of ${latex(system)}$', pad=20)
plt.axhline(y=0, color='black', linestyle='dotted', linewidth=1)
if grid:
plt.grid(True)
if show:
plt.show()
return
return plt
nichols_plot(tf1)