Skip to content

Commit

Permalink
Merge branch 'main' into main
Browse files Browse the repository at this point in the history
  • Loading branch information
@rafapalacios
rafapalacios authored Jun 30, 2023
2 parents df44fce + 328558b commit f4b5472
Show file tree
Hide file tree
Showing 14 changed files with 3,484 additions and 14 deletions.
93 changes: 93 additions & 0 deletions docs/source/content/example_notebooks/UDP_control/control_design_script.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
function [success] = control_design_script(route_directory)
addpath(strcat(route_directory,'/matlab_functions/'));
success = false;
%% Define parameters;
case_name = 'pazy_ROM';
output_folder = strcat(route_directory,'/output/',case_name);
output_folder_linear = strcat(output_folder, '/linear_results/');
complex_matrices = true;

%% Reade from SHARPy generated state space model
[state_space_system, eta_ref] = read_SHARPy_state_space_system(...
strcat(case_name, '.linss.h5'), ...
strcat(output_folder, '/savedata/'), ...
complex_matrices...
);

%% Load simulation settings and store in struct
load(strcat(output_folder_linear, 'simulation_parameters.mat'));
n_nodes = uint8(n_nodes);
rigid_body_motions = false;

input_settings = set_input_parameters(u_inf, ...
state_space_system.Ts, ...
num_modes, ...
num_aero_states, ...
rigid_body_motions, ...
simulation_time, ...
num_control_surfaces, ...
n_nodes, ...
control_input_start, ...
gust_input_start);

%% Remove unused input and outputs from the generated ROM in SHARPy and link deflection and its rate
state_space_system = adjust_state_space_system(state_space_system, input_settings);


%% Get gust input

gust = get_1minuscosine_gust_input(...
gust_length, ...
gust_intensity, ...
state_space_system.Ts, ...
u_inf, ...
simulation_time);
%% Configure Simulink Inputs
model_name = "PID_linear_model";
simIn = Simulink.SimulationInput(model_name);
simIn = setVariable(simIn,'D_gain',0);
simIn = setVariable(simIn,'P_gain',0);
simIn = setVariable(simIn,'I_gain',0);
simIn = setVariable(simIn,'state_space_system',state_space_system);
simIn = setVariable(simIn,'input_settings',input_settings);
simIn = setVariable(simIn,'gust', gust);
%% Run Simulink with PID Controller
warning('off','all');
% open loop
fprintf('Simulate linear open-loop gust response.')
simIn = setVariable(simIn,'controller_on',0);

out_open_loop = sim(simIn);
% closed loop P = 5
simIn = setVariable(simIn,'controller_on',1);
simIn = setVariable(simIn,'P_gain',5);
fprintf('Simulate linear closed-loop gust response with P=5.')
out_PID_5 = sim(simIn);

% closed loop P = 10
simIn = setVariable(simIn,'P_gain',10);
fprintf('Simulate linear closed-loop gust response with P=10.')
out_PID_10 = sim(simIn);

%% Save results
deflection = [out_PID_5.delta.Data(:,1)...
out_PID_10.delta.Data(:,1)...
zeros(size(out_PID_5.delta.Time, 1), 1)];
deflection_rate = [out_PID_5.delta_dot.Data(:,1)...
out_PID_10.delta_dot.Data(:,1)...
zeros(size(out_PID_5.delta_dot.Time, 1), 1)];
tip_deflection = [out_PID_5.actual_output.Data(:,1)...
out_PID_10.actual_output.Data(:,1)...
out_open_loop.actual_output.Data(:,1)];
tip_deflection= tip_deflection + eta_ref(n_nodes/2*6-3);
time_array = out_open_loop.tout;

output_folder_linear = strcat(output_folder, '/linear_results/')
writematrix(deflection, strcat(output_folder_linear, 'deflection.txt'),'Delimiter',',');
writematrix(time_array, strcat(output_folder_linear, 'time_array_linear.txt'),'Delimiter',',');
writematrix(deflection_rate,strcat(output_folder_linear, './deflection_rate.txt'),'Delimiter',',');
writematrix(tip_deflection,strcat(output_folder_linear, './tip_deflection.txt'),'Delimiter',',');

success = true

end
207 changes: 207 additions & 0 deletions docs/source/content/example_notebooks/UDP_control/get_settings_udp.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,207 @@
"""
Script to set the required SHARPy settings for all cases within
the UDP Control example notebook.
"""
def get_settings_udp(model, flow, **kwargs):
gust = kwargs.get('gust', False)

horseshoe = kwargs.get('horseshoe', False)
gravity = kwargs.get('gravity', True)
wake_length_factor = kwargs.get('wake_length_factor', 10)

num_cores = kwargs.get('num_cores', 2)
num_modes = kwargs.get('num_modes', 20)
free_flight = kwargs.get('free_flight', False)
tolerance = kwargs.get('tolerance', 1e-6)
n_load_steps = kwargs.get('n_load_steps', 5)
fsi_tolerance = kwargs.get('fsi_tolerance', 1e-6)
relaxation_factor = kwargs.get('relaxation_factor', 0.1)
newmark_damp = kwargs.get('newmark_damp', 0.5e-4)
output_folder = kwargs.get('output_folder', './output/')
model.config['SHARPy'] = {'case': model.case_name,
'route': model.route,
'flow': flow,
'write_screen': kwargs.get('write_screen', 'on'),
'write_log': 'on',
'save_settings': 'on',
'log_folder': output_folder + '/',
'log_file': model.case_name + '.log'}


model.config['BeamLoader'] = {'unsteady': 'off',
'orientation': model.quat}

model.config['AerogridLoader'] = {'unsteady': 'on',
'aligned_grid': 'on',
'mstar': wake_length_factor*model.M, #int(20/tstep_factor),
'wake_shape_generator': 'StraightWake',
'wake_shape_generator_input': {
'u_inf':model.u_inf,
'u_inf_direction': [1., 0., 0.],
'dt':model.dt,
},
}
if horseshoe:
model.config['AerogridLoader']['mstar'] = 1

model.config['StaticCoupled'] = {
'print_info': 'on',
'max_iter': 200,
'n_load_steps': n_load_steps,
'tolerance': 1e-5,
'relaxation_factor': relaxation_factor,
'aero_solver': 'StaticUvlm',
'aero_solver_settings': {
'rho': model.rho,
'print_info': 'off',
'horseshoe': 'on',
'num_cores': num_cores,
'n_rollup': 0,
'velocity_field_generator': 'SteadyVelocityField',
'velocity_field_input': {
'u_inf': model.u_inf,
'u_inf_direction': model.u_inf_direction},
'vortex_radius': 1e-9},
'structural_solver': 'NonLinearStatic',
'structural_solver_settings': {'print_info': 'off',
'max_iterations': 200,
'num_load_steps': 5,
'delta_curved': 1e-6,
'min_delta': 1e-8,
'gravity': gravity,
'gravity': 9.81},
}


model.config['AerogridPlot'] = {'include_rbm': 'off',
'include_applied_forces': 'on',
'minus_m_star': 0}
model.config['BeamPlot'] = {'include_rbm': 'off',
'include_applied_forces': 'on'}

model.config['WriteVariablesTime'] = {'structure_variables': ['pos'],
'structure_nodes': list(range(0, model.num_node_surf)),
'cleanup_old_solution': 'on',
}



model.config['DynamicCoupled'] = {'print_info': 'on',
'structural_substeps': 10,
'dynamic_relaxation': 'on',
'cleanup_previous_solution': 'on',
'structural_solver': 'NonLinearDynamicPrescribedStep',
'structural_solver_settings': {'print_info': 'off',
'max_iterations': 950,
'delta_curved': 1e-1,
'min_delta': tolerance,
'newmark_damp': newmark_damp,
'gravity': gravity,
'gravity': 9.81,
'num_steps': model.n_tstep,
'dt':model.dt,
},
'aero_solver': 'StepUvlm',
'aero_solver_settings': {'print_info': 'on',
'num_cores': num_cores,
'convection_scheme': 2,
'velocity_field_generator': 'SteadyVelocityField',
'velocity_field_input': {'u_inf': model.u_inf,
'u_inf_direction': [1., 0., 0.]},
'rho': model.rho,
'n_time_steps': model.n_tstep,
'vortex_radius': 1e-9,
'dt': model.dt,
'gamma_dot_filtering': 3},
'fsi_substeps': 200,
'fsi_tolerance': fsi_tolerance,
'relaxation_factor': model.relaxation_factor,
'minimum_steps': 1,
'relaxation_steps': 150,
'final_relaxation_factor': 0.0,
'n_time_steps': model.n_tstep,
'dt': model.dt,
'include_unsteady_force_contribution': kwargs.get('unsteady_force_distribution', True),
'postprocessors': ['WriteVariablesTime', 'BeamPlot', 'AerogridPlot'],
'postprocessors_settings': {'BeamPlot': {'include_rbm': 'on',
'include_applied_forces': 'on'},
'StallCheck': {},
'AerogridPlot': {
'u_inf': model.u_inf,
'include_rbm': 'on',
'include_applied_forces': 'on',
'minus_m_star': 0},
'WriteVariablesTime': {
'structure_variables': ['pos', 'psi'],
'structure_nodes': [model.num_node_surf - 1,
model.num_node_surf,
model.num_node_surf + 1],
},
},
'network_settings': kwargs.get('network_settings', {}),
}

if gust:
gust_settings = kwargs.get('gust_settings', {'gust_shape': '1-cos',
'gust_length': 10.,
'gust_intensity': 0.01,
'gust_offset': 0.})
model.config['DynamicCoupled']['aero_solver_settings']['velocity_field_generator'] = 'GustVelocityField'
model.config['DynamicCoupled']['aero_solver_settings']['velocity_field_input'] = {'u_inf': model.u_inf,
'u_inf_direction': [1., 0, 0],
'relative_motion': bool(not free_flight),
'offset': gust_settings['gust_offset'],
'gust_shape': gust_settings['gust_shape'],
'gust_parameters': {
'gust_length': gust_settings['gust_length'],
'gust_intensity': gust_settings['gust_intensity'] * model.u_inf,
}
}

model.config['PickleData'] = {}

model.config['LinearAssembler'] = {'linear_system': 'LinearAeroelastic',
'inout_coordinates': 'nodes',
# 'recover_accelerations': True,
'linear_system_settings': {
'beam_settings': {'modal_projection': True,
'inout_coords': 'modes',
'discrete_time': True,
'newmark_damp': newmark_damp,
'discr_method': 'newmark',
'dt': model.dt,
'proj_modes': 'undamped',
'num_modes': num_modes,
'print_info': 'on',
'gravity': gravity,
'remove_dofs': []},
'aero_settings': {'dt': model.dt,
'integr_order': 2,
'density': model.rho,
'remove_predictor': True,
'use_sparse': 'off',
'gust_assembler': 'LeadingEdge',
'ScalingDict': kwargs.get('scaling_dict', {'length':1, 'speed': 1, 'density': 1}),
},
'track_body': free_flight,
'use_euler': free_flight,
}}
model.config['Modal'] = {'print_info': True,
'use_undamped_modes': True,
'NumLambda': num_modes,
'rigid_body_modes': free_flight,
'write_modes_vtk': False,
'print_matrices': False,
'continuous_eigenvalues': 'off',
'dt': model.dt,
'plot_eigenvalues': False,
}
model.config['SaveData']['save_linear'] = True
if kwargs.get('remove_gust_input_in_statespace', False):
model.config['LinearAssembler']['linear_system_settings']['aero_settings']['remove_inputs'] = ['u_gust']

rom_settings = kwargs.get('rom_settings', {'use': False})
if rom_settings['use']:
model.config['LinearAssembler']['linear_system_settings']['aero_settings']['rom_method'] = rom_settings['rom_method'],
model.config['LinearAssembler']['linear_system_settings']['aero_settings']['rom_method_settings'] = rom_settings['rom_method_settings']
return model.config
Binary file added docs/source/content/example_notebooks/UDP_control/matlab_functions/PID_linear_model.slx
View file
Binary file not shown.
55 changes: 55 additions & 0 deletions ...source/content/example_notebooks/UDP_control/matlab_functions/adjust_state_space_system.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
function state_space_converted = adjust_state_space_system(state_space_system, input_settings)
%convert_state_space_for_LQR Adjust extracted state space system
% The loaded state space system exported from SHARPy is adjusted here for
% the closed-loop simulations. This includes:
% - extracts the column from the B and D matrices that describe the
% effect of the gust input to the states
% - removing all unused inputs and rearranging for example the
% deflection of a control surface and its rate as both depend on
% each other
% - removes not considered outputs
% - assembles updated discrete state space system


%% Reduce model by deleting unused outupts (only tip deflection relevant here)
C_sensor = state_space_system.C(input_settings.index.tip_displacement,:);
D_sensor = state_space_system.D(input_settings.index.tip_displacement,:);

%% Reduce model with not used inputs
idx_input_end =input_settings.index.control_input_start + 2 * input_settings.num_control_surfaces - 1;

idx_inputs = [input_settings.index.control_input_start:idx_input_end];
B_cs = state_space_system.B(:,idx_inputs);
D_cs = D_sensor(:,idx_inputs);


%% Extract delta to state

new_A_colum = B_cs(:,1 :input_settings.num_control_surfaces);
A = [state_space_system.A new_A_colum; ...
zeros(input_settings.num_control_surfaces,...
(size(state_space_system.A,2)+ input_settings.num_control_surfaces))];

for counter = 0:1:input_settings.num_control_surfaces-1
A(size(A,1)-counter,size(A,2)-counter) = 1; % Explain one
end

%Delete delta input and add delta_dot influence on delta
B_cs(:,1:input_settings.num_control_surfaces) = [];
B_cs = [B_cs; eye(input_settings.num_control_surfaces) * state_space_system.Ts];

% Adding feed through of the delta-input on the output on C and deleting
% column added in C out of D; new column
C = [C_sensor D_cs(:,1:input_settings.num_control_surfaces)];
D_cs(:, 1:input_settings.num_control_surfaces) = [];


%% Get Disturbance Matrices
G = [state_space_system.B(:,input_settings.index.gust_input_start); ...
zeros(input_settings.num_control_surfaces, 1)];
H = D_sensor(:,input_settings.index.gust_input_start);
%% Assemble final state space system

state_space_converted = ss(A,[B_cs G],C,[D_cs H],state_space_system.Ts);

end
Loading

0 comments on commit f4b5472

Please sign in to comment.