hm/sunny360vfe_proj
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
develop
sunny360vfe_proj/tasks/gen_trim_report.m
T
LZY\79495 4570cac9a9 更新子库，固件生成可以
2025-06-04 14:13:22 +08:00

674 lines
14 KiB
Matlab
Raw Permalink Blame History

function gen_trim_report()
%% 计算条件设定
h0 = 50;
gamma0 = 0;
v0 = 10:1:35;
rpt_name = 'trimming_report_20210519';
%% 模型初始化
proj = slproject.getCurrentProject;
wd = cd;
cd([proj.RootFolder '/docs']);
mdl = evalin('base', 'mdl_name');
fcn = [mdl, '_an'];
opts = trimoptions(fcn);
opts.v0_type = uint8(0);
opts.max_alpha = 12/57.3;
opts.min_alpha = -4/57.3;
opts.max_de = 25/57.3;
opts.min_de = -25/57.3;
opts.max_da = 25/57.3;
opts.min_da = -25/57.3;
opts.max_dr = 25/57.3;
opts.min_dr = -25/57.3;
opts.throttle = 0.6;
opts.min_throttle = 0;
opts.TrimOptions.OptimizationOptions.Display = 'iter';
opts.TrimOptions.DisplayReport = 'on';
load_system(fcn);
mdlWks = get_param(fcn,'ModelWorkspace');
%% 初始化报告
import mlreportgen.report.*
import mlreportgen.dom.*
rpt = Report(rpt_name,'docx','temp');
tp = TitlePage;
tp.Title = 'Trimming Report';
tp.Subtitle = ['for model ' fcn];
tp.Author = 'Lyu Weiye';
add(rpt,tp);
%% 计算
trim_lon_report(rpt,opts,mdlWks,h0,v0,gamma0);
% linLon_report(rpt,opts,mdlWks,h0,v0,gamma0);
% trimGamma_report(rpt,opts,mdlWks,h0,v0);
% trimPullup_report(rpt,opts,mdlWks,h0,v0);
% trimSideslip_report(rpt,opts,mdlWks,h0,v0,gamma0);
% trimTurnPhiDr_report(rpt,opts,mdlWks,h0,v0);
%trimTurnPhi_report(rpt,opts,mdlWks,h0,v0);
%% 结束
rptview(rpt);
close_system(fcn, 0); % close without saving
cd(wd);
end
%% 定直平飞配平
function trim_lon_report(rpt,opts,~,h0,v0,gamma0)
hf = figure();
set(gcf,'Position',[0 0 700 600]);
plotTrimLon(hf, ...
opts, ...
h0, v0, ...
gamma0);
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('trimlon(H=%.0f)',h0);
add(rpt,fig);
savefig('trimlon');
delete(gcf);
end
function plotTrimLon(hf, opts, h0, v0, gamma0, xlbl)
if nargin < 6
xlbl = 'V / m/s';
end
j=1;
for v =v0
[success,opts1] = trimLon( opts, h0, v, gamma0 );
if success
opts = opts1;
v_lst(j) = v;
aoa_lst(j) = opts.alpha;
de_lst(j) = opts.de;
thr_lst(j) = opts.throttle;
thrust_lst(j) = opts.thrust;
LoD_lst(j) = opts.LoD;
RPM_lst(j) = opts.RPM;
j=j+1;
end
end
if j>1
figure(hf);
subplot(321);
plot(v_lst, aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^{o}');
hold on;
grid on;
subplot(322);
plot(v_lst, LoD_lst);
xlabel(xlbl);
ylabel('LoD');
hold on;
grid on;
subplot(323);
plot(v_lst, thrust_lst);
xlabel(xlbl);
ylabel('T_r / N');
hold on;
grid on;
subplot(324);
plot(v_lst, de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^{o}');
hold on;
grid on;
subplot(325);
plot(v_lst, RPM_lst);
xlabel(xlbl);
ylabel('RPM');
hold on;
grid on;
end
end
%% 模态分析
function linLon_report(rpt,opts,~,h0,v0,gamma0)
hf = figure();
plotLinLonLat(hf, ...
opts, ...
h0, v0, ...
gamma0);
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('linearize_modes(H=%.0f)',h0);
add(rpt,fig);
savefig('linearize_modes');
delete(gcf);
end
function plotLinLonLat(hf, opts, h0, v0, gamma0)
figure(hf);
for v =v0
[success,opts1,op_point] = trimLon( opts, h0, v, gamma0 );
if success
opts = opts1;
[Glon,Glat] = linfdm(opts.fcn,op_point);
ev = eig(Glon);
subplot(211);
plot(real(ev),imag(ev), 'x');
xlabel('Real');
ylabel('Imag');
legend('lon');
hold on;
ev = eig(Glat);
subplot(212);
plot(real(ev),imag(ev), 'x');
xlabel('Real');
ylabel('Imag');
legend('lat');
hold on;
end
end
end
%% 爬升
function trimGamma_report(rpt,opts,~,h0,v0)
hf = figure();
plotTrimGamma(hf, ...
opts, ...
h0, v0, ...
0);
plotTrimGamma(hf, ...
opts, ...
h0, v0, ...
0.7);
plotTrimGamma(hf, ...
opts, ...
h0, v0, ...
0.8);
plotTrimGamma(hf, ...
opts, ...
h0, v0, ...
0.9);
plotTrimGamma(hf, ...
opts, ...
h0, v0, ...
1);
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('climb_glide(H=%.0f)',h0);
add(rpt,fig);
savefig('climb_glide');
delete(gcf);
end
function plotTrimGamma(hf, opts, h0, v0, thr0, xlbl)
if nargin < 6
xlbl = 'V / m/s';
end
j=1;
for v =v0
[success,opts1] = trimGamma( opts, h0, v, thr0 );
if success && opts1.gamma > -0.2
opts = opts1;
v_lst(j) = v;
aoa_lst(j) = opts.alpha;
de_lst(j) = opts.de;
thr_lst(j) = opts.throttle;
gamma_lst(j) = opts.gamma;
hdot_lst(j) = sin(opts.gamma)*opts.TAS;
thrust_lst(j) = opts.thrust;
RPM_lst(j) = opts.RPM;
j=j+1;
end
end
if j>1
figure(hf);
subplot(321);
plot(v_lst, aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^{o}');
hold on;
grid on;
subplot(322);
plot(v_lst, de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^{o}');
hold on;
grid on;
subplot(323);
plot(v_lst, gamma_lst*57.3);
xlabel(xlbl);
ylabel('\gamma / ^{o}');
hold on;
grid on;
subplot(324);
plot(v_lst, hdot_lst);
xlabel(xlbl);
ylabel('RoC / m/s');
hold on;
grid on;
subplot(325);
plot(v_lst, thrust_lst);
xlabel(xlbl);
ylabel('T / N');
hold on;
grid on;
subplot(326);
plot(v_lst, RPM_lst);
xlabel(xlbl, 'FontSize',12);
ylabel('RPM', 'FontSize',12);
hold on;
grid on;
end
end
%% 定常拉起
function trimPullup_report(rpt,opts,~,h0,v0)
hf = figure();
idx = 1;
for Nz=[-1:0.5:2]
plotTrimPullup(hf, ...
opts, ...
h0, v0, ...
Nz);
idx = idx+1;
end
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('Pullup(H=%.0f)',h0);
add(rpt,fig);
savefig('Pullup');
delete(gcf);
end
function plotTrimPullup(hf, opts, h0, v0, nz0, xlbl)
if nargin < 6
xlbl = 'V / m/s';
end
j=1;
v_lst = [];
for v =v0
[success,opts1] = trimPullup( opts, h0, v, nz0 );
if success
opts = opts1;
v_lst(j) = opts.v;
aoa_lst(j) = opts.alpha;
de_lst(j) = opts.de;
thr_lst(j) = opts.throttle;
nz_lst(j) = opts.Nn;
gamma_lst(j) = opts.gamma;
thrust_lst(j) = opts.thrust;
RPM_lst(j) = opts.RPM;
j=j+1;
end
end
if ~isempty(v_lst)
figure(hf);
subplot(321);
plot(v_lst,aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^o');
grid on
hold on
subplot(322);
plot(v_lst,de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^o');
grid on
hold on
subplot(323);
plot(v_lst, thrust_lst);
xlabel(xlbl);
ylabel('T / N');
grid on
hold on
subplot(324);
plot(v_lst, RPM_lst);
xlabel(xlbl);
ylabel('RPM');
grid on
hold on;
subplot(325);
plot(v_lst, nz_lst);
xlabel(xlbl);
ylabel('Nz');
grid on
hold on;
end
end
%% 抗侧风
function trimSideslip_report(rpt,opts,~,h0,v0,gamma0)
hf = figure();
set(gcf,'Position',[0 0 700 600]);
for beta0 = [ -8 -4 0 4 8 ]
plotTrimSideslip(hf, opts, h0, v0, beta0./57.3, gamma0);
end
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('Sideslip(H=%.0f)',h0);
add(rpt,fig);
savefig('Sideslip');
delete(gcf);
end
function plotTrimSideslip(hf, opts, h0, v0, beta0, gamma0, xlbl)
if nargin <= 6
xlbl = 'V / m/s';
end
j=1;
v_lst = [];
for v =v0
[success,opts1] = trimSideslip( opts, h0, v, beta0, gamma0 );
if success
opts = opts1;
v_lst(j) = opts.v;
aoa_lst(j) = opts.alpha;
aos_lst(j) = opts.beta;
dr_lst(j) = opts.dr;
da_lst(j) = opts.da;
de_lst(j) = opts.de;
phi_lst(j) = opts.phi;
thr_lst(j) = opts.throttle;
thrust_lst(j) = opts.thrust;
RPM_lst(j) = opts.RPM;
j=j+1;
end
end
if ~isempty(v_lst)
figure(hf);
subplot(321);
plot(v_lst, aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^{o}');
hold on;
grid on;
subplot(322);
plot(v_lst, de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^{o}');
hold on;
grid on;
subplot(323);
plot(v_lst, aos_lst*57.3);
xlabel(xlbl);
ylabel('\beta / ^{o}');
hold on;
grid on;
subplot(324);
plot(v_lst, dr_lst*57.3);
xlabel(xlbl);
ylabel('\delta_r / ^{o}');
hold on;
grid on;
subplot(325);
plot(v_lst, da_lst*57.3);
xlabel(xlbl);
ylabel('\delta_a / ^{o}');
hold on;
grid on;
subplot(326);
plot(v_lst, RPM_lst);
xlabel(xlbl);
ylabel('RPM');
hold on;
grid on;
end
end
%% 协调盘旋
function trimTurnPhi_report(rpt,opts,~,h0,v0)
hf = figure();
for phi0 = [-60:15:60]
plotTrimTurnPhi(hf, opts, h0, v0, phi0/57.3);
end
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('turn(H=%.0f)',h0);
add(rpt,fig);
savefig('turn');
delete(gcf);
end
function plotTrimTurnPhi(hf, opts, h0, v0, phi0, xlbl)
if nargin < 6
xlbl = 'V / m/s';
end
j=1;
v_lst = [];
for v = v0
[success,opts1] = trimTurnPhi( opts, h0, v, phi0 );
if success
opts = opts1;
v_lst(j) = opts.v;
aoa_lst(j) = opts.alpha;
aos_lst(j) = opts.beta;
de_lst(j) = opts.de;
da_lst(j) = opts.da;
dr_lst(j) = opts.dr;
thr_lst(j) = opts.throttle;
nz_lst(j) = opts.Nn;
RPM_lst(j) = opts.RPM;
phi_lst(j)= opts.attitude(1);
omega_lst(j) = sqrt(dot(opts.pqr,opts.pqr));
j=j+1;
end
end
if ~isempty(v_lst)
figure(hf);
subplot(331);
plot(v_lst,aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^o');
grid on
hold on
subplot(332);
plot(v_lst,de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^o');
grid on
hold on
subplot(333);
plot(v_lst,RPM_lst);
xlabel(xlbl);
ylabel('RPM');
grid on
hold on
subplot(334);
plot(v_lst,da_lst*57.3);
xlabel(xlbl);
ylabel('\delta_a / ^o');
grid on
hold on
subplot(335);
plot(v_lst,dr_lst*57.3);
xlabel(xlbl);
ylabel('\delta_r / ^o');
grid on
hold on
subplot(336);
plot(v_lst,aos_lst*57.3 );
xlabel(xlbl);
ylabel('\beta / ^o');
grid on
hold on;
subplot(337);
plot(v_lst,phi_lst*57.3);
xlabel(xlbl);
ylabel('\phi / ^o');
grid on
hold on;
subplot(338);
plot(v_lst,nz_lst );
xlabel(xlbl);
ylabel('N_z');
grid on
hold on;
subplot(339);
plot(v_lst,omega_lst*57.3 );
xlabel(xlbl);
ylabel('\oemga / ^o/s');
grid on
hold on;
end
end
%% 非协调盘旋特性
function trimTurnPhiDr_report(rpt,opts,~,h0,v0)
hf = figure();
for phi= -45: 15: 45
plotTurnPhiDr(hf, ...
opts, ...
h0, v0, ...
phi/57.3, 0);
end
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('turn_dr0(H=%.0f)',h0);
add(rpt,fig);
savefig('turn_dr0');
delete(gcf);
end
function plotTurnPhiDr(hf, opts, h0, v0, phi0, dr0, xlbl)
if nargin < 7
xlbl = 'V / m/s';
end
j=1;
for v =v0
[success,opts1] = trimTurnPhiDr( opts, h0, v, phi0, dr0);
if success
opts = opts1;
v_lst(j) = opts.v;
aoa_lst(j) = opts.alpha;
aos_lst(j) = opts.beta;
de_lst(j) = opts.de;
da_lst(j) = opts.da;
dr_lst(j) = opts.dr;
thr_lst(j) = opts.throttle;
nz_lst(j) = opts.Nn;
RPM_lst(j) = opts.RPM;
phi_lst(j)= opts.attitude(1);
omega_lst(j) = sqrt(dot(opts.pqr,opts.pqr));
j=j+1;
end
end
if j>1
figure(hf);
subplot(331);
plot(v_lst,aoa_lst*57.3);
xlabel(xlbl);
ylabel('\alpha / ^o');
grid on
hold on
subplot(332);
plot(v_lst,de_lst*57.3);
xlabel(xlbl);
ylabel('\delta_e / ^o');
grid on
hold on
subplot(333);
plot(v_lst,RPM_lst);
xlabel(xlbl);
ylabel('RPM');
grid on
hold on
subplot(334);
plot(v_lst,da_lst*57.3);
xlabel(xlbl);
ylabel('\delta_a / ^o');
grid on
hold on
subplot(335);
plot(v_lst,dr_lst*57.3);
xlabel(xlbl);
ylabel('\delta_r / ^o');
grid on
hold on
subplot(336);
plot(v_lst,aos_lst*57.3);
xlabel(xlbl);
ylabel('\beta / ^o');
grid on
hold on;
subplot(337);
plot(v_lst,phi_lst*57.3);
xlabel(xlbl);
ylabel('\phi / ^o');
grid on
hold on;
subplot(338);
plot(v_lst,nz_lst);
xlabel(xlbl);
ylabel('N_z');
grid on
hold on;
subplot(339);
plot(v_lst,omega_lst*57.3);
xlabel(xlbl);
ylabel('\omega / ^o/s');
grid on
hold on;
end
end
Reference in New Issue View Git Blame Copy Permalink