diff --git a/ROSCO_toolbox/control_interface.py b/ROSCO_toolbox/control_interface.py index 1a9aceb3..5514bf6a 100644 --- a/ROSCO_toolbox/control_interface.py +++ b/ROSCO_toolbox/control_interface.py @@ -85,7 +85,6 @@ def init_discon(self): # Torque initial condition self.avrSWAP[22] = 0 - # Code this as first call self.avrSWAP[0] = 0 @@ -98,7 +97,6 @@ def init_discon(self): # 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.avcOUTNAME = (self.sim_name + '.RO.dbg').encode('utf-8') self.avcMSG = create_string_buffer(1000) self.discon.DISCON.argtypes = [POINTER(c_float), POINTER(c_int32), c_char_p, c_char_p, c_char_p] # (all defined by ctypes) @@ -106,7 +104,7 @@ def init_discon(self): # Run DISCON self.call_discon() - # Code as not first run + # Code as not first run now that DISCON has been initialized self.avrSWAP[0] = 1 @@ -173,7 +171,6 @@ def call_controller(self, turbine_state, end=False): except KeyError: self.avrSWAP[82] = 0 - # call controller self.call_discon() @@ -182,9 +179,6 @@ def call_controller(self, turbine_state, end=False): self.torque = self.avrSWAP[46] self.nac_yawrate = self.avrSWAP[47] - # if turbine_state['iStatus'] == -1: - # self.avrSWAP[0] = -1 - return(self.torque,self.pitch,self.nac_yawrate) def show_control_values(self): @@ -304,6 +298,9 @@ def __init__(self, network_addresses=["tcp://*:5555", "tcp://*:5556"], verbose=verbose) def get_measurements(self): + ''' + Get measurements from zmq servers + ''' measurements = [None for _ in range(self.nturbs)] for ti in range(self.nturbs): measurements[ti] = self.zmq_servers[ti].get_measurements() @@ -311,7 +308,18 @@ def get_measurements(self): def send_setpoints(self, genTorques=None, nacelleHeadings=None, bladePitchAngles=None): + ''' + Send setpoints to DLL via zmq server for farm level controls + Parameters: + ----------- + genTorques: List + List of generator torques of length self.nturbs + nacelleHeadings: List + List of nacelle headings of length self.nturbs + bladePitchAngles: List + List of blade pitch angles of length self.nturbs + ''' # Default choices if unspecified if genTorques is None: genTorques = [0.0] * self.nturbs @@ -357,6 +365,9 @@ def __init__(self, network_address="tcp://*:5555", identifier="0", self._connect() def _connect(self): + ''' + Connect to zmq server + ''' address = self.network_address # Connect socket @@ -369,11 +380,17 @@ def _connect(self): print("[%s] Successfully established connection with %s" % (self.identifier, address)) def _disconnect(self): + ''' + Disconnect from zmq server + ''' self.socket.close() context = zmq.Context() context.term() def get_measurements(self): + ''' + Receive measurements from ROSCO .dll + ''' if self.verbose: print("[%s] Waiting to receive measurements from ROSCO..." % (self.identifier)) @@ -390,7 +407,6 @@ def get_measurements(self): # Convert to individual strings and then to floats measurements = message_in - # measurements = bytes.decode(message_in) measurements = measurements.replace('\x00', '').split(',') measurements = [float(m) for m in measurements] @@ -421,6 +437,18 @@ def get_measurements(self): def send_setpoints(self, genTorque=0.0, nacelleHeading=0.0, bladePitch=[0.0, 0.0, 0.0]): + ''' + Send setpoints to ROSCO .dll ffor individual turbine control + + Parameters: + ----------- + genTorques: float + Generator torque setpoint + nacelleHeadings: float + Nacelle heading setpoint + bladePitchAngles: List (len=3) + Blade pitch angle setpoint + ''' # Create a message with setpoints to send to ROSCO message_out = b"%016.5f, %016.5f, %016.5f, %016.5f, %016.5f" % ( genTorque, nacelleHeading, bladePitch[0], bladePitch[1], diff --git a/ROSCO_toolbox/controller.py b/ROSCO_toolbox/controller.py index e9ad8ce0..b4e3c353 100644 --- a/ROSCO_toolbox/controller.py +++ b/ROSCO_toolbox/controller.py @@ -10,10 +10,11 @@ # specific language governing permissions and limitations under the License. import numpy as np -import sys, os +import os import datetime -from scipy import interpolate, gradient, integrate +from scipy import interpolate, integrate from ROSCO_toolbox.utilities import list_check +from scipy import optimize # Some useful constants now = datetime.datetime.now() @@ -140,8 +141,6 @@ def __init__(self, controller_params): if self.OL_Mode: ol_params = controller_params['open_loop'] - - self.OL_Ind_Breakpoint = ol_params['OL_Ind_Breakpoint'] self.OL_Ind_BldPitch = ol_params['OL_Ind_BldPitch'] self.OL_Ind_GenTq = ol_params['OL_Ind_GenTq'] @@ -204,15 +203,14 @@ def tune_controller(self, turbine): # separate TSRs by operations regions TSR_below_rated = [min(turbine.TSR_operational, rated_rotor_speed*R/v) for v in v_below_rated] # below rated - TSR_above_rated = rated_rotor_speed*R/v_above_rated # above rated - # TSR_below_rated = np.minimum(np.max(TSR_above_rated), TSR_below_rated) - TSR_op = np.concatenate((TSR_below_rated, TSR_above_rated)) # operational TSRs + TSR_above_rated = rated_rotor_speed*R/v_above_rated # above rated + TSR_op = np.concatenate((TSR_below_rated, TSR_above_rated)) # operational TSRs # Find expected operational Cp values - Cp_above_rated = turbine.Cp.interp_surface(0,TSR_above_rated[0]) # Cp during rated operation (not optimal). Assumes cut-in bld pitch to be 0 + Cp_above_rated = turbine.Cp.interp_surface(0,TSR_above_rated[0]) # Cp during rated operation (not optimal). Assumes cut-in bld pitch to be 0 Cp_op_br = np.ones(len(v_below_rated)) * turbine.Cp.max # below rated Cp_op_ar = Cp_above_rated * (TSR_above_rated/TSR_rated)**3 # above rated - Cp_op = np.concatenate((Cp_op_br, Cp_op_ar)) # operational CPs to linearize around + Cp_op = np.concatenate((Cp_op_br, Cp_op_ar)) # operational CPs to linearize around pitch_initial_rad = turbine.pitch_initial_rad TSR_initial = turbine.TSR_initial @@ -225,22 +223,25 @@ def tune_controller(self, turbine): Ct_op = np.empty(len(TSR_op)) # ------------- Find Linearized State "Matrices" ------------- # + # At each operating point for i in range(len(TSR_op)): - # Find pitch angle as a function of expected operating CP for each TSR + # Find pitch angle as a function of expected operating CP for each TSR operating point Cp_TSR = np.ndarray.flatten(turbine.Cp.interp_surface(turbine.pitch_initial_rad, TSR_op[i])) # all Cp values for a given tsr Cp_maxidx = Cp_TSR.argmax() - Cp_op[i] = np.clip(Cp_op[i], np.min(Cp_TSR[Cp_maxidx:]), np.max(Cp_TSR[Cp_maxidx:])) # saturate Cp values to be on Cp surface # Find maximum Cp value for this TSR - f_cp_pitch = interpolate.interp1d(Cp_TSR[Cp_maxidx:],pitch_initial_rad[Cp_maxidx:]) # interpolate function for Cp(tsr) values - # expected operation blade pitch values + Cp_op[i] = np.clip(Cp_op[i], np.min(Cp_TSR[Cp_maxidx:]), np.max(Cp_TSR[Cp_maxidx:])) # saturate Cp values to be on Cp surface # Find maximum Cp value for this TSR + f_cp_pitch = interpolate.interp1d(Cp_TSR[Cp_maxidx:],pitch_initial_rad[Cp_maxidx:]) # interpolate function for Cp(tsr) values + + # expected operational blade pitch values. Saturates by min_pitch if it exists if v[i] <= turbine.v_rated and isinstance(self.min_pitch, float): # Below rated & defined min_pitch pitch_op[i] = min(self.min_pitch, f_cp_pitch(Cp_op[i])) - elif isinstance(self.min_pitch, float): + elif isinstance(self.min_pitch, float): # above rated & defined min_pitch pitch_op[i] = max(self.min_pitch, f_cp_pitch(Cp_op[i])) - else: + else: # no defined minimum pitch schedule pitch_op[i] = f_cp_pitch(Cp_op[i]) - dCp_beta[i], dCp_TSR[i] = turbine.Cp.interp_gradient(pitch_op[i],TSR_op[i]) # gradients of Cp surface in Beta and TSR directions - dCt_beta[i], dCt_TSR[i] = turbine.Ct.interp_gradient(pitch_op[i],TSR_op[i]) # gradients of Cp surface in Beta and TSR directions + # Calculate Cp Surface gradients + dCp_beta[i], dCp_TSR[i] = turbine.Cp.interp_gradient(pitch_op[i],TSR_op[i]) + dCt_beta[i], dCt_TSR[i] = turbine.Ct.interp_gradient(pitch_op[i],TSR_op[i]) # Thrust Ct_TSR = np.ndarray.flatten(turbine.Ct.interp_surface(turbine.pitch_initial_rad, TSR_op[i])) # all Cp values for a given tsr @@ -259,12 +260,11 @@ def tune_controller(self, turbine): dCt_dbeta = dCt_beta/np.diff(pitch_initial_rad)[0] dCt_dTSR = dCt_TSR/np.diff(TSR_initial)[0] - # Linearized system derivatives - dtau_dbeta = Ng/2*rho*Ar*R*(1/TSR_op)*dCp_dbeta*v**2 - dtau_dlambda = Ng/2*rho*Ar*R*v**2*(1/(TSR_op**2))*(dCp_dTSR*TSR_op - Cp_op) + # Linearized system derivatives, equations from https://wes.copernicus.org/articles/7/53/2022/wes-7-53-2022.pdf + dtau_dbeta = Ng/2*rho*Ar*R*(1/TSR_op)*dCp_dbeta*v**2 # (26) + dtau_dlambda = Ng/2*rho*Ar*R*v**2*(1/(TSR_op**2))*(dCp_dTSR*TSR_op - Cp_op) # (7) dlambda_domega = R/v/Ng dtau_domega = dtau_dlambda*dlambda_domega - dlambda_dv = -(TSR_op/v) Pi_beta = 1/2 * rho * Ar * v**2 * dCt_dbeta @@ -418,23 +418,19 @@ def tune_flap_controller(self,turbine): # Find blade aerodynamic coefficients v_rel = [] phi_vec = [] - alpha=[] for i, _ in enumerate(self.v): turbine.cc_rotor.induction_inflow=True # Axial and tangential inductions try: - a, ap, alpha0, cl, cd = turbine.cc_rotor.distributedAeroLoads( + a, ap, _, _, _ = turbine.cc_rotor.distributedAeroLoads( self.v[i], self.omega_op[i], self.pitch_op[i], 0.0) except ValueError: - loads, derivs = turbine.cc_rotor.distributedAeroLoads( + loads, _ = turbine.cc_rotor.distributedAeroLoads( self.v[i], self.omega_op[i], self.pitch_op[i], 0.0) - a = loads['a'] - ap = loads['ap'] - alpha0 = loads['alpha'] - cl = loads['Cl'] - cd = loads['Cd'] + a = loads['a'] # Axial induction factor + ap = loads['ap'] # Tangential induction factor - # Relative windspeed + # Relative windspeed along blade span v_rel.append([np.sqrt(self.v[i]**2*(1-a)**2 + self.omega_op[i]**2*turbine.span**2*(1-ap)**2)]) # Inflow wind direction phi_vec.append(self.pitch_op[i] + turbine.twist*deg2rad) @@ -449,9 +445,11 @@ def tune_flap_controller(self,turbine): Cdm = np.zeros(num_af) for i,section in enumerate(turbine.af_data): - # assume airfoil section as AOA of zero for slope calculations - for now + # assume airfoil section as AOA of zero for slope calculations a0_ind = section[0]['Alpha'].index(np.min(np.abs(section[0]['Alpha']))) # Coefficients + # - If the flap exists in this blade section, define Cx-plus,-minus,-neutral(0) + # - IF teh flap does not exist in this blade section, Cx matrix is all the same value if section[0]['NumTabs'] == 3: # sections with 3 flaps Clm[i,] = section[0]['Cl'][a0_ind] Cdm[i,] = section[0]['Cd'][a0_ind] @@ -465,12 +463,12 @@ def tune_flap_controller(self,turbine): Cd0[i,] = Cdp[i,] = Cdm[i,] = section[0]['Cd'][a0_ind] Ctrl = float(section[0]['Ctrl']) - # Find slopes + # Find lift and drag coefficient slopes w.r.t. flap angle Kcl = (Clp - Cl0)/( (Ctrl_flp-Ctrl)*deg2rad ) Kcd = (Cdp - Cd0)/( (Ctrl_flp-Ctrl)*deg2rad ) # Find integrated constants - self.kappa = np.zeros(len(v_rel)) + self.kappa = np.zeros(len(v_rel)) # "flap efficacy term" C1 = np.zeros(len(v_rel)) C2 = np.zeros(len(v_rel)) for i, (v_sec,phi) in enumerate(zip(v_rel, phi_vec)): @@ -481,7 +479,6 @@ def tune_flap_controller(self,turbine): # PI Gains if (self.flp_kp_norm == 0 or self.flp_tau == 0) or (not self.flp_kp_norm or not self.flp_tau): raise ValueError('flp_kp_norm and flp_tau must be nonzero for Flp_Mode >= 1') - self.Kp_flap = self.flp_kp_norm / self.kappa self.Ki_flap = self.flp_kp_norm / self.kappa / self.flp_tau @@ -510,22 +507,20 @@ def peak_shaving(self,controller, turbine): ''' # Re-define Turbine Parameters for shorthand - J = turbine.J # Total rotor inertial (kg-m^2) - rho = turbine.rho # Air density (kg/m^3) - R = turbine.rotor_radius # Rotor radius (m) - A = np.pi*R**2 # Rotor area (m^2) - Ng = turbine.Ng # Gearbox ratio (-) - rated_rotor_speed = turbine.rated_rotor_speed # Rated rotor speed (rad/s) + rho = turbine.rho # Air density (kg/m^3) + R = turbine.rotor_radius # Rotor radius (m) + A = np.pi*R**2 # Rotor area (m^2) # Initialize some arrays Ct_op = np.empty(len(controller.TSR_op),dtype='float64') Ct_max = np.empty(len(controller.TSR_op),dtype='float64') - beta_min = np.empty(len(controller.TSR_op),dtype='float64') - # Find unshaved rotor thurst coefficients and associated rotor thrusts - # for i in len(controller.TSR_op): + + # Find unshaved rotor thrust coefficients at each TSR for i in range(len(controller.TSR_op)): Ct_op[i] = turbine.Ct.interp_surface(controller.pitch_op[i],controller.TSR_op[i]) - T = 0.5 * rho * A * controller.v**2 * Ct_op + + # Thrust vs. wind speed + T = 0.5 * rho * A * controller.v**2 * Ct_op # Define minimum max thrust and initialize pitch_min Tmax = controller.ps_percent * np.max(T) @@ -534,18 +529,21 @@ def peak_shaving(self,controller, turbine): # Modify pitch_min if max thrust exceeds limits for i in range(len(controller.TSR_op)): # Find Ct values for operational TSR - # Ct_tsr = turbine.Ct.interp_surface(turbine.pitch_initial_rad, controller.TSR_op[i]) Ct_tsr = turbine.Ct.interp_surface(turbine.pitch_initial_rad,controller.TSR_op[i]) # Define max Ct values Ct_max[i] = Tmax/(0.5 * rho * A * controller.v[i]**2) if T[i] > Tmax: Ct_op[i] = Ct_max[i] else: + # TSR_below_rated = np.minimum(np.max(TSR_above_rated), TSR_below_rated) Ct_max[i] = np.minimum( np.max(Ct_tsr), Ct_max[i]) + # Define minimum pitch angle + # - find min(\beta) so that Ct <= Ct_max and \beta > \beta_fine at each operational TSR f_pitch_min = interpolate.interp1d(Ct_tsr, turbine.pitch_initial_rad, kind='linear', bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1])) pitch_min[i] = max(controller.min_pitch, f_pitch_min(Ct_max[i])) + # Save to controller object controller.ps_min_bld_pitch = pitch_min # save some outputs for analysis or future work @@ -557,9 +555,20 @@ def peak_shaving(self,controller, turbine): self.T = T def min_pitch_saturation(self, controller, turbine): - + ''' + Minimum pitch saturation in low wind speeds to maximize power capture + + Parameters: + ----------- + controller: class + Controller class containing controller operational information + turbine: class + Turbine class containing necessary wind turbine information for controller tuning + ''' # Find TSR associated with minimum rotor speed TSR_at_minspeed = (controller.pc_minspd) * turbine.rotor_radius / controller.v_below_rated + + # For each below rated wind speed operating point for i in range(len(TSR_at_minspeed)): if TSR_at_minspeed[i] > controller.TSR_op[i]: controller.TSR_op[i] = TSR_at_minspeed[i] @@ -567,18 +576,16 @@ def min_pitch_saturation(self, controller, turbine): # Initialize some arrays Cp_op = np.empty(len(turbine.pitch_initial_rad),dtype='float64') min_pitch = np.empty(len(TSR_at_minspeed),dtype='float64') - - # Find Cp-maximizing minimum pitch schedule - # Find Cp coefficients at below-rated tip speed ratios + # ------- Find Cp-maximizing minimum pitch schedule --------- + # Cp coefficients at below-rated tip speed ratios Cp_op = turbine.Cp.interp_surface(turbine.pitch_initial_rad,TSR_at_minspeed[i]) - Cp_max = max(Cp_op) - # f_pitch_min = interpolate.interp1d(Cp_op, -turbine.pitch_initial_rad, kind='quadratic', bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1])) + + # Setup and run small optimization problem to find blade pitch angle that maximizes Cp at a given TSR + # - Finds \beta to satisfy max( Cp(\beta,TSR_op) ) f_pitch_min = interpolate.interp1d(turbine.pitch_initial_rad, -Cp_op, kind='quadratic', bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1])) - from scipy import optimize res = optimize.minimize(f_pitch_min, 0.0) min_pitch[i] = res.x[0] - # min_pitch[i] = f_pitch_min(Cp_max) # modify existing minimum pitch schedule controller.ps_min_bld_pitch[i] = np.maximum(controller.ps_min_bld_pitch[i], min_pitch[i]) @@ -817,19 +824,28 @@ def write_input(self,ol_filename): return open_loop -# helper functions - +# ----------- Helper Functions ----------- def sigma(tt,t0,t1,y0=0,y1=1): ''' generates timeseries for a smooth transition from y0 to y1 from x0 to x1 - inputs: tt - time indices - t0 - start time - t1 - end time - y0 - start output - y1 - end output - - outputs: yy - output timeseries corresponding to tt + Parameters: + ----------- + tt: List-like + time indices + t0: float + start time + t1: float + end time + y0: float + start output + y1: + end output + + Returns: + -------- + yy: List-like + output timeseries corresponding to tt ''' a3 = 2/(t0-t1)**3 @@ -854,14 +870,19 @@ def multi_sigma(xx,x_bp,y_bp): Parameters: ----------- - xx : list of floats (-) + xx: list of floats (-) new sample points - x_bp : list of floats (-) + x_bp: list of floats (-) breakpoints y_bp : list of floats (-) function value at breakpoints + Returns: + -------- + yy: List-like + output timeseries corresponding to tt ''' + # initialize yy = np.zeros_like(xx) # interpolate sigma functions between all breakpoints diff --git a/ROSCO_toolbox/linear/lin_util.py b/ROSCO_toolbox/linear/lin_util.py index 74b2e537..969ea3fa 100644 --- a/ROSCO_toolbox/linear/lin_util.py +++ b/ROSCO_toolbox/linear/lin_util.py @@ -102,21 +102,33 @@ def pc_sensitivity(linturb, controller, u): return sens def smargin(linturb, controller, u_eval): - # try: - # k_float = inputs[1] - # except: + ''' + Calculates the stability margin for an open-loop system + + linturb: object + LinearTurbineModel object + controller: object + ROSCO toolbox controller object + u_eval: float + wind speed to evaluate system at + ''' + + # Find standard transfer functions sens_sys = pc_sensitivity(linturb, controller, u_eval) ol_sys = pc_openloop(linturb, controller, u_eval) - cl_sys = pc_closedloop(linturb, controller, u_eval) + # Convert to state space sp_plant = sp.signal.StateSpace(ol_sys.A, ol_sys.B, ol_sys.C, ol_sys.D) sp_sens = sp.signal.StateSpace(sens_sys.A, sens_sys.B, sens_sys.C, sens_sys.D) + # Minimum distance to the critical point on the nyquist diagram def nyquist_min(om): return np.abs(sp.signal.freqresp(sp_plant, w=om)[1] + 1.) + # Maximum value of sensitivity function def sens_min(om): return -sp.signal.bode(sp_sens, w=om)[1] with warnings.catch_warnings(): - warnings.simplefilter("ignore") - # Find first local maxima in sensitivity function + warnings.simplefilter("ignore") # NJA: Lots of scipy errors because of poorly posed state matrices - ignore them + + # Find first local maximum in bode magnitude of the sensitivity function ws, m, _ = sp.signal.bode(sp_sens, n=100000) m0 = m[0] m1 = m[1] @@ -128,35 +140,45 @@ def sens_min(om): return -sp.signal.bode(sp_sens, w=om)[1] i += 1 w0 = ws[i] + # --- Find important magnitudes and related frequencies with just the sampled values --- + # # NJA - optimization is run in the next "step" to fine tune these values further + # magnitude of sensitivity function at first local maximum sm_mag = sp.signal.freqresp(sp_plant, w=w0)[1] sm = np.sqrt((1 - np.abs(sm_mag.real))**2 + sm_mag.imag**2) + # magnitude and frequency of distance to critical point on nyquist nearest_nyquist = nyquist_min(ws).min() nearest_nyquist_freq = ws[nyquist_min(ws).argmin()] mag_at_min = sp.signal.freqresp(sp_plant, w=nearest_nyquist_freq)[1] + + # If any poles of the sensitivity function are unstable, solve an optimization problem on the Nyquist trajectory + # to minimize the distance to the critical point over all frequencies. Start the gradient-based optimization at + # the nearest value calculated in the sweep above + #: NJA - optimization is done to reduce precision errors if any(sp_sens.poles > 0): if nearest_nyquist < sm: res = sp.optimize.minimize(nyquist_min, nearest_nyquist_freq, method='SLSQP', options={ 'finite_diff_rel_step': 1e-8}) # Make sure this didn't fail if res.status != 0: + # if optimization failed, just use the closest value from the initial sweep sm2 = nearest_nyquist else: sm2 = min(abs(res.fun), abs(nearest_nyquist)) + # Dubious error catching because of strange nyquist shapes in unstable systems sm_list = [sm, sm2] mag_list = [np.abs(sm_mag), np.abs(mag_at_min)] sm = sm_list[np.argmax(mag_list)] - sm *= -1 # Flip sign because it's unstable - else: + sm *= -1 # Flip sign to have a "negative sensitivity margin" because it's unstable + else: # system is stable res = sp.optimize.minimize(nyquist_min, nearest_nyquist_freq, method='SLSQP', options={'finite_diff_rel_step': 1e-6}) if res.status != 0: + # if optimization failed, just use the closest value from the initial sweep sm = nearest_nyquist else: sm = min(res.fun, nearest_nyquist) - - return sm def interp_plant(linturb, v, return_scipy=True): @@ -247,6 +269,7 @@ def add_pcomp(linturb, k_float): Modified linturb with parallel compensation term included ''' + # Find relevant state and input indices state_str = 'derivative of 1st tower fore-aft' state_idx = np.flatnonzero(np.core.defchararray.find( linturb.DescStates, state_str) > -1).tolist() @@ -254,8 +277,10 @@ def add_pcomp(linturb, k_float): input_idx = np.flatnonzero(np.core.defchararray.find( linturb.DescCntrlInpt, input_str) > -1).tolist() + # Modify linear system to be \dot{x} = (A+B_fK)x + Bu, + # - B_fKx accounts for parallel compensation in the linear system K = np.zeros((linturb.B_ops.shape[1], linturb.A_ops.shape[1], linturb.A_ops.shape[2])) - K[input_idx, state_idx, :] = -k_float # NJA: negative to account of OF linearization sign conventions + K[input_idx, state_idx, :] = -k_float # NJA: negative to account for OpenFAST linearization sign conventions linturb2 = copy.copy(linturb) linturb2.A_ops = linturb.A_ops + linturb.B_ops * K diff --git a/ROSCO_toolbox/linear/robust_scheduling.py b/ROSCO_toolbox/linear/robust_scheduling.py index bb347bab..70bc418e 100644 --- a/ROSCO_toolbox/linear/robust_scheduling.py +++ b/ROSCO_toolbox/linear/robust_scheduling.py @@ -1,5 +1,6 @@ ''' Methods for finding robust gain schedules +- Implemented in the OpenMDAO framework ''' import numpy as np @@ -38,7 +39,7 @@ def setup(self): 'Error: omega_pc and zeta_pc must be scalars for robust controller tuning.') # Load ROSCO Turbine and Controller - if 'dict_inputs' in ROSCO_options.keys(): # Allow for turbine parameters to be passed in as a dictionary + if 'dict_inputs' in ROSCO_options.keys(): # Allow for turbine parameters to be passed in as a dictionary (from WEIS) dict_inputs = ROSCO_options['dict_inputs'] # Define turbine based on inputs self.turbine = type('', (), {})() @@ -79,7 +80,7 @@ def setup(self): self.controller = ROSCO_controller.Controller(ROSCO_options['controller_params']) self.controller.tune_controller(self.turbine) - else: + else: # otherwise define controller and turbine objection self.controller, self.turbine = load_ROSCO(ROSCO_options['path_params'], ROSCO_options['turbine_params'], ROSCO_options['controller_params']) @@ -101,6 +102,9 @@ def setup(self): desc='Maximized controller bandwidth') def compute(self, inputs, outputs): + ''' + Computes the stability margin for a given controller bandwidth (omega) + ''' k_float = inputs['k_float'][0] if k_float: linturb = add_pcomp(self.linturb, k_float) @@ -121,10 +125,14 @@ def compute(self, inputs, outputs): class rsched_driver(): ''' - A driver for scheduling robust controllers + A driver for scheduling robust controllers. + Wrapper for the RobustScheduling OpenMDAO Explicit Component + + - Example 12 shows how to use this ''' def __init__(self, options): + # initialize options self.linturb_options = options['linturb_options'] self.ROSCO_options = options['ROSCO_options'] self.path_options = options['path_options'] @@ -136,6 +144,7 @@ def __init__(self, options): self.output_dir = self.path_options['output_dir'] self.output_name = self.path_options['output_name'] + # setup levels, especially useful for DOEs if 'levels' in self.opt_options.keys(): self.opt_options['levels'] = self.opt_options['levels'] else: @@ -162,9 +171,12 @@ def setup(self): else: self.om_problem = om.Problem() + # Add robust scheduling subsystem self.om_problem.model.add_subsystem('r_sched', RobustScheduling(linturb_options=self.linturb_options, ROSCO_options=self.ROSCO_options)) + # Setup design of experiments or optimization problem. + # - NJA: DOEs can only be done for one wind speed at a time if self.opt_options['driver'] == 'design_of_experiments': self.om_problem = self.init_doe(self.om_problem, levels=self.opt_options['levels']) if list_check(self.opt_options['windspeed']): @@ -199,12 +211,8 @@ def setup(self): self.om_doe = self.om_problem self.om_doe = self.add_dv(self.om_doe, ['omega', 'k_float']) if self.opt_options['driver'] == 'optimization': - # self.om_doe = copy.deepcopy(self.om_problem) - # self.om_doe = self.add_dv(self.om_doe, ['omega']) - self.om_opt = self.om_problem self.om_opt = self.add_dv(self.om_opt, ['omega', 'k_float']) - # self.om_opt = self.init_optimization(self.om_opt) def execute(self): ''' @@ -222,10 +230,12 @@ def execute(self): self.omegas = [] self.k_floats = [] self.sms = [] + + # Iterate optimization over each wind speed in opt_options for u in self.opt_options['windspeed']: - om0 = 0.1 + om0 = 0.1 # Initial omega - hard coded try_count = 0 - while try_count < 3: + while try_count < 3: # Allow this to try three times before failing. NJA: sometimes small initial condition changes can improve convergence # Setup optimization self.om_opt.set_val('r_sched.u_eval', u) self.om_opt.set_val('r_sched.omega', om0) @@ -240,10 +250,11 @@ def execute(self): self.om_opt.cleanup() if self.om_opt.driver.fail: + # Restart with a new initial omega if the optimizer failed try_count += 1 om0 = np.random.random_sample(1)*self.opt_options['omega'][-1] else: - try_count = 4 + try_count += 1 # save values self.omegas.append(self.om_opt.get_val('r_sched.omega')[0]) @@ -251,6 +262,9 @@ def execute(self): self.sms.append(self.om_opt.get_val('r_sched.sm')[0]) def plot_schedule(self): + ''' + Plots tuning value and stability margins w.r.t. wind speed + ''' fig, ax = plt.subplots(3, 1, constrained_layout=True, sharex=True) ax[0].plot(self.opt_options['windspeed'], self.omegas) ax[0].set_ylabel('omega_pc') @@ -265,7 +279,17 @@ def plot_schedule(self): plt.show() def add_dv(self, om_problem, opt_vars): - '''add design variables''' + ''' + Add design variables to OM problem + + Parameters: + ----------- + om_problem: om.problem + RobustScheduling OpenMDAO problem + opt_vars: list + Variables to optimize, e.g.:['omega', 'k_float'] + + ''' if 'omega' in opt_vars and list_check(self.opt_options['omega']): om_problem.model.add_design_var( @@ -276,6 +300,7 @@ def add_dv(self, om_problem, opt_vars): 'r_sched.k_float', lower=self.opt_options['k_float'][0], upper=self.opt_options['k_float'][1], ref=100) # Make sure design variables are stored appropriately in OM problem + # - NJA: This is mostly needed for WEIS integration as a nested OpenMDAO problem om_problem.model._design_vars = om_problem.model._static_design_vars return om_problem @@ -297,8 +322,6 @@ def init_doe(self, om_problem, levels=20): # om_problem.driver = om.DOEDriver(om.LatinHypercubeGenerator(samples=levels)) # om_problem.driver = om.DOEDriver(om.UniformGenerator(num_samples=20)) os.makedirs(self.output_dir, exist_ok=True) - # om_problem.driver.options['run_parallel'] = True - # om_problem.driver.options['procs_per_model'] = 1 return om_problem @staticmethod @@ -334,6 +357,23 @@ def post_doe(self, save_csv=False): def load_DOE(doe_logs, outfile_name=None): + ''' + Loads and processes doe log files + - OpenMDAO DOEs generate a large set of log files, this function collects them + + Parameters: + ----------- + doe_logs: str,list + string (single) or list (multiple) of doe log file names + outfile_name: str, optional + name of output .csv file to save data + + Returns: + -------- + df: pd.DataFrame + Pandas dataframe containing collected DOE data + + ''' if isinstance(doe_logs, str): doe_logs = [doe_logs] @@ -376,6 +416,7 @@ def load_DOE(doe_logs, outfile_name=None): def load_OMsql(log): + '''load OpenMDAO sql file''' print('loading {}'.format(log)) cr = om.CaseReader(log) rec_data = {} @@ -391,6 +432,16 @@ def load_OMsql(log): def load_linturb(linfile_path, load_parallel=False): + ''' + Load linear turbine models + + Parameters: + ----------- + linfile_path: string + Path to folder containing .lin files + load_parllel: bool + Load parallel True/False + ''' # Parse openfast linearization filenames filenames = glob.glob(os.path.join(linfile_path, '*.lin')) linfiles = [os.path.split(file)[1] for file in filenames] @@ -404,6 +455,18 @@ def load_linturb(linfile_path, load_parallel=False): def load_ROSCO(path_params, turbine_params, controller_params): + ''' + Load ROSCO controller and turbine objects + + Parameters: + ----------- + path_params: dict + Path parameters from tuning yaml + turbine_params: dict + Turbine parameters from tuning yaml + controller_params: dict + Controller parameters from tuning yaml + ''' turbine = ROSCO_turbine.Turbine(turbine_params) controller = ROSCO_controller.Controller(controller_params) @@ -415,7 +478,6 @@ def load_ROSCO(path_params, turbine_params, controller_params): if __name__ == '__main__': - # Setup linear turbine paths linfile_path = os.path.join(os.path.dirname(os.path.dirname( os.path.dirname(os.path.abspath(__file__)))), 'Test_Cases', 'IEA-15-240-RWT-UMaineSemi', 'linearizations') diff --git a/ROSCO_toolbox/sim.py b/ROSCO_toolbox/sim.py index 6fcfdb86..54e3df62 100644 --- a/ROSCO_toolbox/sim.py +++ b/ROSCO_toolbox/sim.py @@ -106,9 +106,6 @@ def sim_ws_series(self, t_array, ws_array, rotor_rpm_init=10, init_pitch=0.0, ws = ws_array[i] wd = wd_array[i] - # ws_x = ws_array[i]*np.cos(wd_array[i] - yaw_angle) - # ws_y = ws_array[i]*np.sin(wd_array[i] - yaw_angle) - # Load current Cq data tsr = rot_speed[i-1] * self.turbine.rotor_radius / ws cq = self.turbine.Cq.interp_surface(bld_pitch[i-1], tsr) @@ -116,15 +113,12 @@ def sim_ws_series(self, t_array, ws_array, rotor_rpm_init=10, init_pitch=0.0, # 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**3) * (cp/tsr) * ws**2 - # 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.turbine.GenEff/100 - self.turbine.Ng * gen_torque[i-1]) gen_speed[i] = rot_speed[i] * self.turbine.Ng # -- Simple nacelle model nac_yawerr[i] = wd - nac_yaw[i-1] - # print('rot_speed = {}'.format(cp)) - # populate turbine state dictionary turbine_state = {} if i < len(t_array): @@ -142,7 +136,7 @@ def sim_ws_series(self, t_array, ws_array, rotor_rpm_init=10, init_pitch=0.0, turbine_state['Yaw_fromNorth'] = nac_yaw[i] turbine_state['Y_MeasErr'] = nac_yawerr[i-1] - + # Define outputs gen_torque[i], bld_pitch[i], nac_yawrate[i] = self.controller_int.call_controller(turbine_state) # Calculate the power diff --git a/ROSCO_toolbox/turbine.py b/ROSCO_toolbox/turbine.py index 0176a663..e97d7203 100644 --- a/ROSCO_toolbox/turbine.py +++ b/ROSCO_toolbox/turbine.py @@ -151,11 +151,13 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_ txt_filename: str, optional filename for *.txt, only used if rot_source='txt' """ + # Use weis if it exists if use_weis: from weis.aeroelasticse.FAST_reader import InputReader_OpenFAST else: from ROSCO_toolbox.ofTools.fast_io.FAST_reader import InputReader_OpenFAST + # Load OpenFAST model using the FAST_reader print('Loading FAST model: %s ' % FAST_InputFile) self.TurbineName = FAST_InputFile.strip('.fst') fast = self.fast = InputReader_OpenFAST() @@ -163,6 +165,7 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_ fast.FAST_directory = FAST_directory fast.execute() + # Use Performance tables if defined, otherwise use defaults if txt_filename: self.rotor_performance_filename = txt_filename else: @@ -175,7 +178,7 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_ self.hubHt = fast.fst_vt['ElastoDyn']['TowerHt'] + fast.fst_vt['ElastoDyn']['Twr2Shft'] self.NumBl = fast.fst_vt['ElastoDyn']['NumBl'] self.TowerHt = fast.fst_vt['ElastoDyn']['TowerHt'] - self.shearExp = 0.2 #HARD CODED FOR NOW + self.shearExp = 0.2 #NOTE: HARD CODED if 'default' in str(fast.fst_vt['AeroDyn15']['AirDens']): fast.fst_vt['AeroDyn15']['AirDens'] = 1.225 self.rho = fast.fst_vt['AeroDyn15']['AirDens'] @@ -196,7 +199,6 @@ def load_from_fast(self, FAST_InputFile,FAST_directory, FAST_ver='OpenFAST',dev_ self.J = self.rotor_inertia + self.generator_inertia * self.Ng**2 self.rated_torque = self.rated_power/(self.GenEff/100*self.rated_rotor_speed*self.Ng) self.rotor_radius = self.TipRad - # self.omega_dt = np.sqrt(self.DTTorSpr/self.J) # Load blade information self.load_blade_info() @@ -267,18 +269,14 @@ def load_from_ccblade(self): pitch_flat = pitch_mesh.flatten() omega_mesh, _ = np.meshgrid(omega_array, pitch_initial) omega_flat = omega_mesh.flatten() - # tsr_flat = (omega_flat * rpm2RadSec * self.rotor_radius) / ws_flat # Get values from cc-blade print('Running CCBlade aerodynamic analysis, this may take a minute...') - try: # wisde/master as of Nov 9, 2020 - _, _, _, _, CP, CT, CQ, CM = self.cc_rotor.evaluate(ws_flat, omega_flat, pitch_flat, coefficients=True) - except(ValueError): # wisdem/dev as of Nov 9, 2020 - outputs, derivs = self.cc_rotor.evaluate(ws_flat, omega_flat, pitch_flat, coefficients=True) - CP = outputs['CP'] - CT = outputs['CT'] - CQ = outputs['CQ'] + outputs, derivs = self.cc_rotor.evaluate(ws_flat, omega_flat, pitch_flat, coefficients=True) + CP = outputs['CP'] + CT = outputs['CT'] + CQ = outputs['CQ'] print('CCBlade aerodynamic analysis run successfully.') # Reshape Cp, Ct and Cq @@ -321,8 +319,7 @@ def generate_rotperf_fast(self, openfast_path, FAST_runDirectory=None, run_BeamD from ROSCO_toolbox.ofTools.case_gen import runFAST_pywrapper, CaseGen_General from ROSCO_toolbox.ofTools.util import FileTools # Load pCrunch tools - from pCrunch import pdTools, Processing - + from pCrunch import Processing # setup values for surface v0 = self.v_rated + 2 @@ -595,14 +592,15 @@ def __init__(self,performance_table, pitch_initial_rad, TSR_initial): if len(self.max_ind[1]) > 1: print('ROSCO_toolbox Warning: repeated maximum values in a performance table and the last one @ pitch = {} rad. was taken...'.format(self.pitch_opt[-1])) - TSR_ind = np.arange(0,len(TSR_initial)) - TSR_fine_ind = np.linspace(TSR_initial[0],TSR_initial[-1],int(TSR_initial[-1] - TSR_initial[0])*100) - f_TSR = interpolate.interp1d(TSR_initial,TSR_initial,bounds_error='False',kind='quadratic') # interpolate function for Cp(tsr) values - TSR_fine = f_TSR(TSR_fine_ind) - f_performance = interpolate.interp1d(TSR_initial,performance_beta_max,bounds_error='False',kind='quadratic') # interpolate function for Cp(tsr) values - performance_fine = f_performance(TSR_fine_ind) - performance_max_ind = np.where(performance_fine == np.max(performance_fine)) - self.TSR_opt = float(TSR_fine[performance_max_ind[0]]) + # Find TSR that maximizes Cx at fine pitch + # - TSR to satisfy: max( Cx(TSR, \beta_fine) ) = TSR_opt + TSR_fine_ind = np.linspace(TSR_initial[0],TSR_initial[-1],int(TSR_initial[-1] - TSR_initial[0])*100) # Range of TSRs to interpolate accross + f_TSR = interpolate.interp1d(TSR_initial,TSR_initial,bounds_error='False',kind='quadratic') # interpolate function for Cp(tsr) values + TSR_fine = f_TSR(TSR_fine_ind) # TSRs at fine pitch + f_performance = interpolate.interp1d(TSR_initial,performance_beta_max,bounds_error='False',kind='quadratic') # interpolate function for Cx(tsr) values + performance_fine = f_performance(TSR_fine_ind) # Cx values at fine pitch + performance_max_ind = np.where(performance_fine == np.max(performance_fine)) # Find max performance at fine pitch + self.TSR_opt = float(TSR_fine[performance_max_ind[0]]) # TSR to maximize Cx at fine pitch def interp_surface(self,pitch,TSR): ''' diff --git a/docs/source/api_change.rst b/docs/source/api_change.rst index 795d0bea..eec47b4b 100644 --- a/docs/source/api_change.rst +++ b/docs/source/api_change.rst @@ -75,23 +75,18 @@ Line Input Name Example Value ROSCO v2.4.1 to ROSCO v2.5.0 ------------------------------- Two filter parameters were added to - - - change the high pass filter in the floating feedback module - - change the low pass filter of the wind speed estimator signal that is used in torque control Open loop control inputs, users must specify: - - The open loop input filename, an example can be found in Examples/Example_OL_Input.dat - - Indices (columns) of values specified in OL_Filename IPC - Proportional Control capabilities were added, 1P and 2P gains should be specified ====== ================= ====================================================================================================================================================================================================== -Line Input Name Example Value +Line Flag Name Example Value ====== ================= ====================================================================================================================================================================================================== 20 OL_Mode 0 ! OL_Mode - Open loop control mode {0: no open loop control, 1: open loop control vs. time, 2: open loop control vs. wind speed} 27 F_WECornerFreq 0.20944 ! F_WECornerFreq - Corner frequency (-3dB point) in the first order low pass filter for the wind speed estimate [rad/s].