Nyquist plotting (#157)

* change default interp type

* improve nested MPI handling

* fix cl bug, sm calc more robust, allow single plant for interp

* add load from pickle capabilites

* add re-try logic

* save v_above_rated to controller object

* fix BAR cp_ct_cq surface path

* update discons

* remove rogue DISCON

* Add linear model visualization tools

* plot nyquist

* remove plt.show()

* fix indexing, windspeed

Examples/example_12.py

@@ -17,7 +17,8 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from ROSCO_toolbox.inputs.validation import load_rosco_yaml
-from ROSCO_toolbox.linear.robust_scheduling import rsched_driver
+from ROSCO_toolbox.linear.robust_scheduling import rsched_driver, load_linturb
+from ROSCO_toolbox.linear.lin_vis import lin_plotting
 from ROSCO_toolbox import turbine as ROSCO_turbine
 from ROSCO_toolbox import controller as ROSCO_controller
 
@@ -123,6 +124,28 @@ def run_example():
     if False:
         plt.show()
     else:
-        plt.savefig(os.path.join(example_out_dir, '12_RobustSched.png'))
+        fig.savefig(os.path.join(example_out_dir, '12_RobustSched.png'))
+
+    # ---- Plot nyquist ----
+    # Re-load and trimlinturb for plotting
+    linturb = load_linturb(options['linturb_options']['linfile_path'], load_parallel=True) # 
+    linturb.trim_system(desInputs=['collective'], desOutputs=['RtSpeed'])
+
+    # Plotting parameters
+    u = 12
+    omega = 0.1
+    k_float = 0.0
+    controller.zeta_pc = controller.zeta_pc[0]
+
+    # plot
+    lv = lin_plotting(controller, turbine, linturb)
+    xlim=ylim=[-2,2]
+    lv.plot_nyquist(u, omega, k_float=k_float, xlim=xlim, ylim=ylim)
+
+    if False:
+        plt.show()
+    else:
+        plt.savefig(os.path.join(example_out_dir, '12_Nyquist.png'))
+
 if __name__ == '__main__':
     run_example()

ROSCO_toolbox/linear/lin_vis.py

@@ -0,0 +1,87 @@
+'''
+Visualization helpers for linear models
+'''
+from ROSCO_toolbox.linear.lin_util import add_pcomp, pc_openloop
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as sp
+
+
+class lin_plotting():
+    def __init__(self, controller, turbine, linturb):
+        '''
+        Parameters
+        ----------
+        controller: object
+            ROSCO controller object
+        turbine: object
+            ROSCO turbine object
+        linturb: object
+            ROSCO linturb object
+        '''
+        self.turbine    = turbine
+        self.controller = controller
+        self.linturb    = linturb
+
+    def plot_nyquist(self, u, omega, k_float=0.0, xlim=None, ylim=None, fig=None, ax=None, num='Nyquist', **kwargs):
+        '''
+        Plot nyquist diagram
+        
+        Parameters:
+        -----------
+        u: float
+            windspeed for linear model
+        omega: float
+            controller bandwidth
+        k_float: float, optional
+            parallel compensation feedback gain
+        '''
+        if fig and ax:
+            self.fig = fig
+            self.ax = ax
+        else:
+            self.fig, self.ax = plt.subplots(1,1,num=num)
+        w, H = self.get_nyquistdata(u, omega, k_float=k_float)
+        self.line, = self.ax.plot(H.real, H.imag,  **kwargs)
+        plt.scatter(-1,0,marker='x',color='r')
+
+        self.ax.set_xlim(xlim)
+        self.ax.set_ylim(ylim)
+        plt.xlabel('Real Axis')
+        plt.ylabel('Imaginary Axis')
+        plt.title('Nyquist Diagram\n$u = {}, \omega = {}, k_f = {}$'.format(u,omega,k_float))
+        plt.grid(True)
+
+
+    def get_nyquistdata(self, u, omega, k_float=0.0):
+        '''
+        Gen nyquist diagram data to plot
+        
+        Parameters:
+        -----------
+        u: float
+            windspeed for linear model
+        omega: float
+            controller bandwidth
+        k_float: float, optional
+            parallel compensation feedback gain
+        '''
+        self.controller.omega_pc = omega
+        self.controller.U_pc = u
+        self.controller.tune_controller(self.turbine)
+        if k_float:
+            linturb = add_pcomp(self.linturb, k_float)
+        else:
+            linturb = self.linturb
+        sys_ol = pc_openloop(linturb, self.controller, u)
+        sys_sp = sp.signal.StateSpace(sys_ol.A, sys_ol.B, sys_ol.C, sys_ol.D)
+        w, H = sp.signal.freqresp(sys_sp)
+
+        return w, H
+
+def sm_circle(r):
+    theta = np.linspace(0,2*np.pi,100)
+    x = r * np.cos(theta) - 1
+    y = r * np.sin(theta)
+
+    return x,y