diff --git a/WTC_toolbox/control_interface.py b/WTC_toolbox/control_interface.py index 626518f87..a729a0a3b 100644 --- a/WTC_toolbox/control_interface.py +++ b/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) \ No newline at end of file + print('Torque',self.torque) \ No newline at end of file diff --git a/WTC_toolbox/sim.py b/WTC_toolbox/sim.py index 20a020210..2e6d58b43 100644 --- a/WTC_toolbox/sim.py +++ b/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() + \ No newline at end of file diff --git a/WTC_toolbox/turbine.py b/WTC_toolbox/turbine.py index 3b18f348d..1b9941952 100644 --- a/WTC_toolbox/turbine.py +++ b/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 diff --git a/examples/example_05.py b/examples/example_05.py index 9a1109f95..d3d6fc7a1 100644 --- a/examples/example_05.py +++ b/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 diff --git a/examples/example_06.py b/examples/example_06.py new file mode 100644 index 000000000..19f72c83a --- /dev/null +++ b/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)