update

explore_wave.py

@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 17 17:42:06 2024
+
+@author: Yihan Liu
+"""
+
+def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases):
+    """
+    This function generates the Pierson-Moskowitz spectrum for a given wind speed U10 and frequency f.
+    
+    parameters
+    ----------
+    U19_5 : float
+        the average wind speed at 19.5m above the sea surface
+    zeta : float
+        the x component to evaluate
+    eta : float
+        the y component to evaluate. (Note: the coordinate system here is different
+                                      from the Betti model. The downward is negative
+                                      in this case)
+    t: float
+        the time to evaluate.
+    random_phase : Numpy Array
+        the random phase to generate wave. Should be in [0, 2*pi)
+
+    Returns
+    -------
+    wave_eta : float
+        The wave elevation
+    [v_x, v_y, a_x, a_y]: list
+        The wave velocity and acceleration in x and y direction
+    """
+    g = 9.81  # gravitational constant
+    alpha = 0.0081  # Phillips' constant
+
+    f_pm = 0.14*(g/U19_5)  # peak frequency
+
+    N = 400
+
+    cutof_f = 3*f_pm # Cutoff frequency
+
+    f = np.linspace(0.1, cutof_f, N) # Array
+    omega = 2*np.pi*f # Array
+    delta_f = f[1] - f[0] # Array
+
+    S_pm = (alpha*g**2/((2*np.pi)**4*f**5))*np.exp(-(5/4)*(f_pm/f)**4) # Array
+
+    a = np.sqrt(2*S_pm*delta_f)
+    k = omega**2/g    
+
+    # Generate random phases all at once
+
+
+    # Perform the calculations in a vectorized manner
+    sin_component = np.sin(omega*t - k*zeta + random_phases)
+    cos_component = np.cos(omega*t - k*zeta + random_phases)
+    exp_component = np.exp(k*eta)
+
+    wave_eta = np.sum(a * sin_component)
+
+    v_x = np.sum(omega * a * exp_component * sin_component)
+    v_y = np.sum(omega * a * exp_component * cos_component)
+
+    a_x = np.sum((omega**2) * a * exp_component * cos_component)
+    a_y = -np.sum((omega**2) * a * exp_component * sin_component)
+
+    return wave_eta, [v_x, v_y, a_x, a_y]

gen_torque_controller.py

@@ -263,8 +263,8 @@ def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases):
     a_x = np.sum((omega**2) * a * exp_component * cos_component)
     a_y = -np.sum((omega**2) * a * exp_component * sin_component)
 
-    #return wave_eta, [v_x, v_y, a_x, a_y]
-    return 0, [0,0,0,0]
+    return wave_eta, [v_x, v_y, a_x, a_y]
+    #return 0, [0,0,0,0]
 
 
 def structure(x_1, beta, omega_R, t, Cp_type, performance, v_w, v_aveg, random_phases):
@@ -572,6 +572,9 @@ def WindTurbine(omega_R, v_in, beta, T_E, t, Cp):
     rho = 1.225 # (kg/m^3) Density of air
     A = 12469 # (m^2) Rotor area
     eta_G = 97 # (-) Speed ratio between high and low speed shafts
+    P_0 = 5296610 # rated power
+
+    T_E = P_0 / (eta_G*omega_R)
 
     tildeJ_R = eta_G**2*J_G + J_R
     tildeT_E = eta_G*T_E
@@ -583,6 +586,8 @@ def WindTurbine(omega_R, v_in, beta, T_E, t, Cp):
     T_A = P_A/omega_R
     domega_R = (1/tildeJ_R)*(T_A - tildeT_E)
 
+    #print(T_E*eta_G*omega_R)
+
     return domega_R
 
 
@@ -666,6 +671,7 @@ def rk4(Betti, x0, t0, tf, dt, beta_0, T_E, Cp_type, performance, v_w, v_wind, s
     """
 
     d_BS = 37.550 # (m) The position of center of weight of BS (platform and tower)
+    eta_G = 97 # (-) Speed ratio between high and low speed shafts
 
     n = int((tf - t0) / dt) + 1
     t = np.linspace(t0, tf, n)
@@ -686,11 +692,12 @@ def rk4(Betti, x0, t0, tf, dt, beta_0, T_E, Cp_type, performance, v_w, v_wind, s
 
     def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i):
 
-
-        eta_G = 97 # (-) Speed ratio between high and low speed shafts
         J_G = 534.116 # (kg*m^2) Total inertia of electric generator and high speed shaft
         J_R = 35444067 # (kg*m^2) Total inertia of blades, hub and low speed shaft
         tildeJ_R = eta_G**2*J_G + J_R
+
+
+
 
         rated_omega_R = 1.26711 # The rated rotor speed is 12.1 rpm
         #rated_omega_R = 1.571
@@ -704,6 +711,9 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i):
         K_p = 0.0765*(2*tildeJ_R*rated_omega_R*zeta_phi*omega_phin*GK)/(eta_G*(-dpdbeta_0))
         K_i = 0.013*(tildeJ_R*rated_omega_R*omega_phin**2*GK)/(eta_G*(-dpdbeta_0))
         K_d = 0.187437
+        #K_p = (2*tildeJ_R*rated_omega_R*zeta_phi*omega_phin*GK)/(eta_G*(-dpdbeta_0))
+        #K_i = (tildeJ_R*rated_omega_R*omega_phin**2*GK)/(eta_G*(-dpdbeta_0))
+        #K_d = 0
 
         error_omega_R = omega_R - rated_omega_R
         error[i] = error_omega_R
@@ -738,6 +748,7 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i):
     h_waves = []
     for i in range(n - 1):
         betas.append(beta)
+
         k1, h_wave, Q_t = Betti(x[i], t[i], beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases)
         k2 = Betti(x[i] + 0.5 * dt * k1, t[i] + 0.5 * dt, beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases)[0]
         k3 = Betti(x[i] + 0.5 * dt * k2, t[i] + 0.5 * dt, beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases)[0]
@@ -1169,7 +1180,7 @@ def plot_fft(wave_eta, t):
 
 
 
-#reproduce_save_driver(seeds)
+reproduce_save_driver(seeds)
 
 
 

reproduce_results.py

@@ -12,6 +12,7 @@
 import matplotlib.pyplot as plt
 from scipy.signal import spectrogram
 from matplotlib.lines import Line2D
+from gen_wind_Van_Der_Hoven import generate_wind
 
 def process_rotor_performance(input_file = "Cp_Ct.NREL5MW.txt"):
     """
@@ -100,6 +101,64 @@ def CpCtCq(TSR, beta, performance):
     # Get the C_p value at the index 
     return C_p[TSR_index][pitch_index], C_t[TSR_index][pitch_index]
 
+def gen_turbulence(v_bar, L, k_sigma_v, T_s, N_t, white_noise, 
+                   delta_omega = 0.002, M = 5000, N = 100):
+    """
+    Generate turbulencec component for wind speed
+
+    Parameters
+    ----------
+    v_bar : int
+        Average wind speed
+    L : int
+        Turbulence length
+    k_sigma_v : float
+        Slope parameter
+    T_s : int
+        Time step
+    N_t : int
+        Number of step
+    white_noise : np.array
+        the white noise with mean = 0 and std = 1, has length 
+
+    Returns
+    -------
+    Array of wind speed with turbulence
+
+    """
+
+    # Step 1: Update the current values of the parameters in 
+    # the turbulence component model
+    T_F = L / v_bar
+    sigma_v = k_sigma_v * v_bar
+    K_F = np.sqrt((2 * np.pi * T_F)/(4.20654 * T_s))
+
+    # Step 2: Calculate the discrete impulse response of the filter
+    delta_omega = 0.002 # Frequency step size
+    M = 5000 # Number of frequency points
+    N = 100 # Numerical parameters for convolution integration, divide 
+            # finite integral from 0 to t to N regions
+
+    # Discrete frequency domain P(omega)
+    P = np.zeros(M + 1)
+    for r in range(M + 1):
+        P[r] = np.real(K_F / (1 + 1j * r * delta_omega * T_F)**(5/6))
+
+    # Discrete impulse response h(k) === h(T_s*k), k range from 0 to N
+    h = np.zeros(N + 1)
+    for k in range(N + 1):
+        h[k] = T_s * delta_omega * (2/np.pi) * np.sum(P * np.cos(k * np.arange(M + 1) * T_s * delta_omega))
+
+    # Step 3: Generate the turbulence component in the interval using convolution
+    v_t = np.zeros(N_t + 1)
+
+    # Zero-pad the white noise 
+    white_noise_padded = np.pad(white_noise, (0, N), 'constant')
+
+    for m in range(N_t + 1):
+        v_t[m] = T_s * np.sum(h * white_noise_padded[m : m + N + 1])
+
+    return v_bar + sigma_v * v_t
 
 def genWind(v_w, end_time, time_step, seed):
     """
@@ -263,8 +322,8 @@ def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases):
     a_x = np.sum((omega**2) * a * exp_component * cos_component)
     a_y = -np.sum((omega**2) * a * exp_component * sin_component)
 
-    #return wave_eta, [v_x, v_y, a_x, a_y]
-    return 0, [0,0,0,0]
+    return wave_eta, [v_x, v_y, a_x, a_y]
+    #return 0, [0,0,0,0]
 
 
 def structure(x_1, beta, omega_R, t, Cp_type, performance, v_w, v_aveg, random_phases):
@@ -772,18 +831,18 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i):
     # dicard data for first 500s
     discard_steps = int(500 / 0.5)
 
-    np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_{seeds[2]}_wind_wave_time.npz', 
+    np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_wind_wave_time.npz', 
                                                     t=t,  
                                                     v_wind=v_wind, 
                                                     wave_eta=wave_eta)    
 
-    t_sub = t[::steps][discard_steps:]
-    x_sub = x[::steps][discard_steps:]
-    v_wind_sub = v_wind[:len(t)][::steps][discard_steps:]
-    wave_eta_sub = np.array(wave_eta)[::steps][discard_steps:]
-    h_wave_sub = np.array(h_waves)[::steps][discard_steps:]
-    betas_sub = betas[::steps][discard_steps:]
-    Qt_list_sub = Qt_list[::steps][discard_steps:]
+    t_sub = t[::steps][:-discard_steps]
+    x_sub = x[::steps][:-discard_steps]
+    v_wind_sub = v_wind[:len(t)][::steps][:-discard_steps]
+    wave_eta_sub = np.array(wave_eta)[::steps][:-discard_steps]
+    h_wave_sub = np.array(h_waves)[::steps][:-discard_steps]
+    betas_sub = betas[::steps][:-discard_steps]
+    Qt_list_sub = Qt_list[::steps][:-discard_steps]
 
     return t_sub-t_sub[0], x_sub, v_wind_sub, wave_eta_sub, h_wave_sub, betas_sub, Qt_list_sub
 
@@ -818,9 +877,25 @@ def main(end_time, v_w, x0, seeds, seed_wave, time_step = 0.05, Cp_type = 0):
     # modify this to change initial condition
     #[zeta, v_zeta, eta, v_eta, alpha, omega, omega_R]
     #v_wind = genWind(v_w, end_time, time_step, seeds)
-    v_wind = np.load(f'reproduced_results/turbsim_output/{seeds[0]}_{seeds[1]}.npy')
+    #v_wind = np.load(f'reproduced_results/turbsim_output/{seeds[0]}_{seeds[1]}.npy')
     #v_wind = genWind_seeds(seeds)
     #v_wind = np.full(45000, 20)
+    def insert_elements(arr, num_elements):
+
+        result = []
+        for i in range(len(arr) - 1):
+            start, end = arr[i], arr[i + 1]
+            result.extend(np.linspace(start, end, num_elements + 2)[:-1])
+        result.append(arr[-1])
+        return np.array(result)
+
+    # generate a random seed
+    state_before = np.random.get_state()
+    np.random.seed(seeds)
+    white_noise = np.random.randn(end_time + 1)
+    np.random.set_state(state_before)
+
+    v_wind = insert_elements(gen_turbulence(v_w, 180, 0.13, 1, end_time, white_noise), int(1/time_step))
 
     # modify this to change run time and step size
     #[Betti, x0 (initial condition), start time, end time, time step, beta, T_E]
@@ -882,8 +957,8 @@ def reproduce_save_driver(seeds):
     v_w = 20
     end_time = 2000 #end_time < 3000
 
-    seeds_wind = [seeds[0], seeds[1]]
-    seed_wave = seeds[2]
+    seeds_wind = seeds[0]
+    seed_wave = seeds[1]
 
 
     x0 = np.array([-2.61426271, 
@@ -896,7 +971,7 @@ def reproduce_save_driver(seeds):
     t, x, v_wind, wave_eta, h_wave, betas, Q_t = main(end_time, v_w, x0, seeds_wind, seed_wave)
     end_time -= 500
 
-    np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_{seeds[2]}.npz', 
+    np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}.npz', 
                                                     t=t,  
                                                     x=x, 
                                                     v_wind=v_wind, 
@@ -916,25 +991,25 @@ def load_data(seeds):
     load the percentile and extreme value
     '''
 
-    output_file_name = f'{seeds[0]}_{seeds[1]}_{seeds[2]}.npz'
+    output_file_name = f'{seeds[0]}_{seeds[1]}.npz'
 
     data = np.load(f'reproduced_results/data/{output_file_name}', allow_pickle=True)
 
     # Extracting the data
     t = data['t'][:-1]
     state = data['x'][:-1]
-    #beta = np.rad2deg(data['betas'])
+    beta = np.rad2deg(data['betas'])
     x = data['x'][:-1]
     wind_speed = data['v_wind'][:-1]
     wave_eta = data['wave_eta'][:-1]
     data.close()
-
+    '''
     pitch_rate = x[:, 5]  
     pitch_acceleration = np.diff(pitch_rate)
     last_acceleration = pitch_acceleration[-1][None]
     pitch_acceleration = np.concatenate((pitch_acceleration, last_acceleration), axis=0)[:, None] 
     state = np.concatenate((x[:, :6], pitch_acceleration), axis=1)
-
+    '''
 
 
     # Extracting percentile data
@@ -985,8 +1060,8 @@ def plot_helper(ax):
         ax[1].set_xlim(0, t[-1])
 
         # plot 7 states
-        for j in range(7):
-        #for j in range(6):
+        #for j in range(7):
+        for j in range(6):
             ax[j+2].plot(t, max_state[:,j], alpha=0.6, color='green', linewidth=0.5)
             ax[j+2].plot(t, min_state[:,j], alpha=0.6, color='orange', linewidth=0.5)
 
@@ -1004,7 +1079,7 @@ def plot_helper(ax):
 
             ax[j+2].tick_params(axis='both', labelsize=16) 
 
-        '''
+
         ax[8].plot(t, state[:, -1], color='black', linewidth=0.5)
         ax[8].set_xlabel('Time (s)', fontsize=12)
         #ax[j+2].set_ylabel(f'{state_names[j]}')
@@ -1037,7 +1112,7 @@ def plot_helper(ax):
                            Line2D([0], [0], color='orange', lw=1, alpha=0.6, label='The Minimum at Each Time Step')]
         
         ax[9].legend(handles=legend_elements, loc='center', fontsize=17.5)
-
+        '''
 
     # for 8 states including pitch acceleration:
 
@@ -1112,7 +1187,7 @@ def plot_fft(wave_eta, t):
 #seeds = [966870, 8017384, 8986784]
 #surge events
 #seeds = [-4966241,  2503014,  9071735]
-seeds = [ -402337, -6699134,  7762480]
+#seeds = [ -402337, -6699134,  7762480]
 #surge velocity
 #2eeds = [-3919707,  4205247,   211083]
 #wave short corr
@@ -1154,7 +1229,7 @@ def plot_fft(wave_eta, t):
 #seeds = [ 4200242, -3440907,  5773012] #std:  0.31226751594476654
 
 #std = load_data(seeds)
-
+seeds =[7619098, 6844915]
 
 seed = [[-8615404,  1149694,  9191470],
         [3973823, 4556159, 3377501],
@@ -1169,7 +1244,7 @@ def plot_fft(wave_eta, t):
 
 
 
-#reproduce_save_driver(seeds)
+reproduce_save_driver(seeds)