Merge pull request #5 from NREL/start_sim

Start sim

WTC_toolbox/control_interface.py

@@ -13,6 +13,11 @@
 import numpy as np
 from numpy.ctypeslib import ndpointer
 
+# Some useful constants
+deg2rad = np.deg2rad(1)
+rad2deg = np.rad2deg(1)
+rpm2RadSec = 2.0*(np.pi)/60.0
+
 class ConInt():
     """
     Define interface to a given controller
@@ -21,6 +26,11 @@ class ConInt():
     def __init__(self, lib_name):
         """
         Setup the interface
+
+        Parameters:
+        -----------
+            lib_name = str
+                        name of dynamic library containing controller, (.dll,.so,.dylib)
         """
         self.lib_name = lib_name
         self.param_name = 'DISCON.IN' # this appears to be hard-coded so match here for now
@@ -35,12 +45,11 @@ def __init__(self, lib_name):
         # Initialize
         self.pitch = 0
         self.torque = 0
-
+        # -- discon
         self.discon = cdll.LoadLibrary(self.lib_name)
-
         self.avrSWAP = np.zeros(self.avr_size)
 
-        # Define avrSWAP some
+        # Define some avrSWAP parameters
         self.avrSWAP[2] = self.DT
         self.avrSWAP[60] = self.num_blade
 
@@ -52,48 +61,56 @@ def __init__(self, lib_name):
         self.avrSWAP[49] = self.char_buffer
         self.avrSWAP[50] = self.char_buffer
 
+        # Initialize DISCON and related
         self.aviFAIL = c_int32() # 1
         self.accINFILE = self.param_name.encode('utf-8')
         self.avcOUTNAME = create_string_buffer(1000) # 'DEMO'.encode('utf-8')
         self.avcMSG = create_string_buffer(1000)
-
-        # self.discon.DISCON.argtypes = [ndpointer(c_double, flags="C_CONTIGUOUS"), POINTER(c_int32), c_char_p, c_char_p, c_char_p] # (all defined by ctypes)
         self.discon.DISCON.argtypes = [POINTER(c_float), POINTER(c_int32), c_char_p, c_char_p, c_char_p] # (all defined by ctypes)
 
+        # Run DISCON
         self.call_discon()
 
-
-        # First call
-        # self.discon.DISCON(self.avrSWAP, byref(self.aviFAIL), self.accINFILE, self.avcOUTNAME, self.avcMSG)
-
-
         # Code as not first run
         self.avrSWAP[0] = 1
 
 
     def call_discon(self):
-
-        # Convert AVR swap to the right thing
+        '''
+        Call DISCON.dll (or .so,.dylib)
+        '''
+        # Convert AVR swap to the c pointer
         c_float_p = POINTER(c_float)
         data = self.avrSWAP.astype(np.float32)
-        data_p = data.ctypes.data_as(c_float_p)
+        p_data = data.ctypes.data_as(c_float_p)
 
-        self.discon.DISCON(data_p, byref(self.aviFAIL), self.accINFILE, self.avcOUTNAME, self.avcMSG)
+        # Run DISCON
+        self.discon.DISCON(p_data, byref(self.aviFAIL), self.accINFILE, self.avcOUTNAME, self.avcMSG)
 
         # Push back to avr swap
-        print(data_p[47])
-        print(self.avrSWAP[47])
-        #for i in range(len(self.avrSWAP)):
-        #    self.avrSWAP
-        #print('len',len(data_p),len(self.avrSWAP))
-
-
+        self.avrSWAP = data
 
 
     def call_controller(self,t,dt,pitch,genspeed,rotspeed,ws):
+        '''
+        Runs the controller. Passes current turbine state to the controller, and returns control inputs back
         
-
-
+        Parameters:
+        -----------
+            t: float
+                time, (s)
+            dt: float
+                timestep, (s)
+            pitch: float
+                blade pitch, (rad)
+            genspeed: float
+                generator speed, (rad/s)
+            rotspeed: float
+                rotor speed, (rad/s)
+            ws: float
+                wind speed, (m/s)
+        '''
+
         # Add states to avr
         self.avrSWAP[1] = t
         self.avrSWAP[2] = dt
@@ -103,15 +120,12 @@ def call_controller(self,t,dt,pitch,genspeed,rotspeed,ws):
         self.avrSWAP[26] = ws
 
         self.call_discon()
-        # self.discon.DISCON(self.avrSWAP, byref(self.aviFAIL), self.accINFILE, self.avcOUTNAME, self.avcMSG)
 
         self.pitch = self.avrSWAP[41]
-        self.torque = self.avrSWAP[47]
+        self.torque = self.avrSWAP[46]
+
+        return(self.torque,self.pitch)
 
     def show_control_values(self):
         print('Pitch',self.pitch)
-        print('Torque',self.torque)
-
-
-#discon.DISCON.argtypes = [POINTER(c_float), c_int, c_char_p, c_char_p, c_char_p] # (all defined by ctypes)
-# discon.DISCON(byref(avrSWAP), aviFAIL, accINFILE, avcOUTNAME, avcMSG)
+        print('Torque',self.torque)

WTC_toolbox/sim.py

@@ -8,3 +8,132 @@
 # under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
 # CONDITIONS OF ANY KIND, either express or implied. See the License for the
 # specific language governing permissions and limitations under the License.
+
+import numpy as np
+from WTC_toolbox import turbine as wtc_turbine
+from WTC_toolbox import control_interface as ci
+import matplotlib.pyplot as plt
+import sys
+
+# Some useful constants
+deg2rad = np.deg2rad(1)
+rad2deg = np.rad2deg(1)
+rpm2RadSec = 2.0*(np.pi)/60.0
+
+class Sim():
+    """
+    Define interface to a given controller
+    """
+
+    def __init__(self, turbine, controller_int):
+        """
+        Setup the simulator
+        """
+        self.turbine = turbine
+        self.controller_int = controller_int
+
+        # TOTAL TEMPORARY HACK
+        self.gb_eff = 0.95
+        self.gen_eff = 0.95
+
+
+    def sim_ws_series(self,t_array,ws_array,rotor_rpm_init=10,init_pitch=0.0, make_plots=True):
+        '''
+        Simulate simplified turbine model using a complied controller (.dll or similar).
+            - currently a 1DOF rotor model
+
+        Parameters:
+        -----------
+            t_array: float
+                     Array of time steps, (s)
+            ws_array: float
+                      Array of wind speeds, (s)
+            rotor_rpm_init: float, optional
+                            initial rotor speed, (rpm)
+            init_pitch: float, optional
+                        initial blade pitch angle, (deg)
+            make_plots: bool, optional
+                        True: generate plots, False: don't. 
+        '''
+        # Store turbine data for conveniente
+        dt = t_array[1] - t_array[0]
+        R = self.turbine.RotorRad
+        GBRatio = self.turbine.Ng
+
+        # Declare output arrays
+        pitch = np.ones_like(t_array) * init_pitch 
+        rot_speed = np.ones_like(t_array) * rotor_rpm_init * rpm2RadSec # represent rot speed in rad / s
+        gen_speed = np.ones_like(t_array) * rotor_rpm_init * GBRatio # represent gen speed in rad/s
+        aero_torque = np.ones_like(t_array) * 1000.0
+        gen_torque = np.ones_like(t_array) # * trq_cont(turbine_dict, gen_speed[0])
+        gen_power = np.ones_like(t_array) * 0.0
+
+        # Test for cc_blade Cq information. 
+        #       - If not, assume available matrices loaded from text file, and interpolate those
+        try: 
+            self.turbine.cc_rotor.evaluate(ws_array[1], [rot_speed[0]/rpm2RadSec], 
+                                                        [pitch[0]], 
+                                                        coefficients=True)
+            use_interpolated = False
+        except: 
+            use_interpolated = True
+            print('CCBLADE Error -- attempting to use interpolated Cp, Ct, Cq, tables')
+
+        # Loop through time
+        for i, t in enumerate(t_array):
+            if i == 0:
+                continue # Skip the first run
+            ws = ws_array[i]
+
+            # Load current Cq data
+            if use_interpolated:
+                tsr = rot_speed[i-1] * self.turbine.RotorRad / ws
+                cq = self.turbine.Cq.interp_surface([pitch[i-1]],tsr)
+            if not use_interpolated:
+                P, T, Q, M, Cp, Ct, cq, CM = self.turbine.cc_rotor.evaluate([ws], 
+                                                        [rot_speed[i-1]/rpm2RadSec], 
+                                                        [pitch[i-1]], 
+                                                        coefficients=True)
+
+            # Update the turbine state
+            #       -- 1DOF model: rotor speed and generator speed (scaled by Ng)
+            aero_torque[i] = 0.5 * self.turbine.rho * (np.pi * R**2) * cq * R * ws**2
+            rot_speed[i] = rot_speed[i-1] + (dt/self.turbine.J)*(aero_torque[i] * self.gb_eff - self.turbine.Ng * gen_torque[i-1])
+            gen_speed[i] = rot_speed[i] * self.turbine.Ng 
+
+            # Call the controller
+            gen_torque[i], pitch[i] = self.controller_int.call_controller(t,dt,pitch[i-1],gen_speed[i],rot_speed[i],ws)
+
+            # Calculate the power
+            gen_power[i] = gen_speed[i] * np.pi/30.0 * gen_torque[i] * self.gen_eff
+
+        # Save these values
+        self.pitch = pitch
+        self.rot_speed = rot_speed
+        self.gen_speed = gen_speed
+        self.aero_torque = aero_torque
+        self.gen_torque = gen_torque
+        self.gen_power = gen_power
+        self.t_array = t_array
+        self.ws_array = ws_array
+
+        if make_plots:
+            fig, axarr = plt.subplots(4,1,sharex=True,figsize=(6,10))
+
+            ax = axarr[0]
+            ax.plot(self.t_array,self.ws_array,label='Wind Speed')
+            ax.grid()
+            ax.legend()
+            ax = axarr[1]
+            ax.plot(self.t_array,self.rot_speed,label='Rot Speed')
+            ax.grid()
+            ax.legend()
+            ax = axarr[2]
+            ax.plot(self.t_array,self.gen_torque,label='Gen Torque')
+            ax.grid()
+            ax.legend()
+            ax = axarr[3]
+            ax.plot(self.t_array,self.pitch,label='Bld Pitch')
+            ax.grid()
+            ax.legend()
+

WTC_toolbox/turbine.py

@@ -124,7 +124,7 @@ def load_from_fast(self, FAST_InputFile,FAST_directory,drivetrain_inertia, FAST_
         self.Ng = fast.fst_vt['ElastoDyn']['GBRatio']
         self.rho = rho
         self.RotorRad = TipRad
-
+        
         # Calculate rated rotor speed for now by scaling from NREL 5MW
         self.RRspeed = (63. / TipRad) * 12.1 * rpm2RadSec
 

examples/example_05.py

@@ -1,5 +1,6 @@
 # Example_05
 # make some test controller calls
+# This example is not so useful and can perhaps be deleted when simulator example is working
 
 
 from WTC_toolbox import turbine as wtc_turbine

examples/example_06.py

@@ -0,0 +1,48 @@
+# Example_06
+# Step wind simulation
+
+from WTC_toolbox import turbine as wtc_turbine
+from WTC_toolbox import sim as wtc_sim
+from WTC_toolbox import control_interface as ci
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+# ensure proper directory location 
+os.chdir('/Users/nabbas/Documents/WindEnergyToolbox/WTC_toolbox/examples')
+
+# Load turbine model
+# Initialiize a turbine class
+turbine = wtc_turbine.Turbine()
+
+# Load turbine from OpenFAST and *.txt file
+FAST_InputFile = '5MW_Land.fst'
+FAST_directory = '/Users/nabbas/Documents/TurbineModels/NREL_5MW/5MW_Land'
+txt_filename = 'Cp_Ct_Cq.txt'
+drivetrain_inertia = 40469564.444
+turbine.load_from_fast(FAST_InputFile,FAST_directory,drivetrain_inertia,dev_branch=True,rot_source='txt', txt_filename=txt_filename)
+
+# # Load turbine quick from python
+# turbine.load('saved_turbine.p')
+
+# Load controller 
+lib_name = 'test_controller/DISCON.dll'
+# lib_name = '/Users/pfleming/Desktop/git_tools/floating/DRC_Fortran/DISCON//DISCON_glin64.so'
+controller_int = ci.ConInt(lib_name)
+
+# Load the simulator
+sim = wtc_sim.Sim(turbine,controller_int)
+
+# Define a wind speed history
+dt = 0.1
+tlen = 400      # length of time to simulate (s)
+ws0 = 9         # initial wind speed (m/s)
+t= np.arange(0,tlen,dt) 
+ws = np.ones_like(t) * ws0
+# add steps at every 100s
+for i in range(len(t)):
+    ws[i] = ws[i] + t[i]//100
+
+
+# Run simulator
+sim.sim_ws_series(t,ws)