#+TITLE: SVD Control :DRAWER: #+STARTUP: overview #+LANGUAGE: en #+EMAIL: dehaeze.thomas@gmail.com #+AUTHOR: Dehaeze Thomas #+HTML_LINK_HOME: ../index.html #+HTML_LINK_UP: ../index.html #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_MATHJAX: align: center tagside: right font: TeX #+PROPERTY: header-args:matlab :session *MATLAB* #+PROPERTY: header-args:matlab+ :comments org #+PROPERTY: header-args:matlab+ :results none #+PROPERTY: header-args:matlab+ :exports both #+PROPERTY: header-args:matlab+ :eval no-export #+PROPERTY: header-args:matlab+ :output-dir figs #+PROPERTY: header-args:matlab+ :tangle no #+PROPERTY: header-args:matlab+ :mkdirp yes #+PROPERTY: header-args:shell :eval no-export #+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}") #+PROPERTY: header-args:latex+ :imagemagick t :fit yes #+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150 #+PROPERTY: header-args:latex+ :imoutoptions -quality 100 #+PROPERTY: header-args:latex+ :results raw replace :buffer no #+PROPERTY: header-args:latex+ :eval no-export #+PROPERTY: header-args:latex+ :exports both #+PROPERTY: header-args:latex+ :mkdirp yes #+PROPERTY: header-args:latex+ :output-dir figs :END: * Simscape Model - Gravimeter ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> #+end_src #+begin_src matlab :exports none :results silent :noweb yes <> #+end_src ** Simulink #+begin_src matlab open('gravimeter.slx') #+end_src #+begin_src matlab %% Name of the Simulink File mdl = 'gravimeter'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/F1'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/F2'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/F3'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Acc_side'], 1, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Acc_side'], 2, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Acc_top'], 1, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Acc_top'], 2, 'openoutput'); io_i = io_i + 1; G = linearize(mdl, io); G.InputName = {'F1', 'F2', 'F3'}; G.OutputName = {'Ax1', 'Az1', 'Ax2', 'Az2'}; #+end_src The plant as 6 states as expected (2 translations + 1 rotation) #+begin_src matlab :results output replace size(G) #+end_src #+RESULTS: : State-space model with 4 outputs, 3 inputs, and 6 states. #+begin_src matlab :exports none freqs = logspace(-2, 2, 1000); figure; for in_i = 1:3 for out_i = 1:4 subplot(4, 3, 3*(out_i-1)+in_i); plot(freqs, abs(squeeze(freqresp(G(out_i,in_i), freqs, 'Hz'))), '-'); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); end end #+end_src #+begin_src matlab :tangle no :exports results :results file replace exportFig('figs/open_loop_tf.pdf', 'width', 'full', 'height', 'full'); #+end_src #+name: fig:open_loop_tf #+caption: Open Loop Transfer Function from 3 Actuators to 4 Accelerometers #+RESULTS: [[file:figs/open_loop_tf.png]] ** Matlab Code :noexport: #+begin_src matlab clc; % close all g = 100000; w0 = 2*pi*.5; % MinusK BM1 tablle l = 0.5; %[m] la = 1; %[m] h = 1.7; %[m] ha = 1.7;% %[m] m = 400; %[kg] k = 15e3;%[N/m] kv = k; kh = 15e3; I = 115;%[kg m^2] % c = 0.06; % l = 0.4719; %[m] % la = .477; %[m] % h = 1.8973; %[m] % ha = 1.9060;% %[m] % m = 98.1421; %[kg] % k = 14512;%[N/m] % I = 28.5372;%[kg m^2] cv = 0.03; ch = 0.03; %% System definition [Fr, x1, z1, x2, z2, wx, wz, x12, z12, PHIwx, PHIwz,xsum,zsum,xdelta,zdelta,rot]... = modelGeneration(m,I,k,h,ha,l,la,cv,ch,kv,kh); %% Bode options P = bodeoptions; P.FreqUnits = 'Hz'; P.MagUnits = 'abs'; P.MagScale = 'log'; P.Grid = 'on'; P.PhaseWrapping = 'on'; P.Xlim = [1e-1,1e2]; %P.PhaseVisible = 'off'; w = 2*pi*logspace(-1,2,1000); %% curves points % slide 4 F_sl4 = [2e-1 4e-1 7e-1 1 2 3 5]; Amp_sl4 = [ 1 2 4 2.5 1 7e-1 7e-1]; F_sl4_phase = [2e-1 4e-1 7e-1 1 ]; Phase_sl4 = (180/pi).*[0 0 -0.5 -1.7]; %slide 6 F_sl6 = [2e-1 4e-1 1 2 3 5]; Amp_sl6 = [1 1 6e-1 2e-1 3e-1 3e-1]; F_sl6_phase = [2e-1 4e-1 1 ]; Phase_sl6 = (180/pi).*[0 0 0 ]; %slide 9 F_sl9 = [2.5e-1 4e-1 6e-1 1 1.7 2.2 3 5 10]; Amp_sl9 = [3 6 1 5e-1 1 2 7e-1 2.5e-1 7e-2]; Phase_sl9 = (180/pi)*[0 -1 -pi 0 -1 -1.5 -pi -pi -pi]; % slide 14 F_sl14 = [ 2e-1 4e-1 6e-1 8e-1 1 2 3 5 10]; Amp_sl14 = [9e-1 1.5 1.2 0.35 .3 1.2 .3 .1 5e-2]; F_sl14_phase = [ 2e-1 4e-1 6e-1 8e-1 ]; Phase_sl14 = (180/pi).*[0 0 -1.7 -2]; %rotation F_rot = [1e-1 2e-1 4e-1 5e-1 7e-1 1 2 3 6.5 10 20]; Amp_rot = [7e-8 2.2e-7 3e-7 1e-7 2e-8 9e-9 3e-8 9e-9 1e-9 4e-10 8e-11]; %% Plots % %slide 3 % figure % loglog(Fr,abs(x2).^.5,Fr,abs(x1).^.5,Fr,abs(xsum).^.5,Fr,abs(xdelta).^.5) % xlabel('Frequency [Hz]');ylabel('Acceleration [m/s^2/rtHz]') % legend('Top sensor','Bottom sensor','Half sum','Half difference'); % title('Horizontal') % xlim([7e-2 1e1]); %slide 4 figure subplot 211 loglog(Fr, abs(x12)./abs(x1));hold on; loglog(F_sl4,Amp_sl4,'*'); xlabel('Frequency [Hz]');ylabel('Amplitude [-]'); title('X direction Top/bottom sensor'); xlim([7e-2 1e1]); subplot 212 semilogx(Fr, (180/pi).*angle(x12./abs(x1)));hold on; loglog(F_sl4_phase,Phase_sl4,'*'); xlabel('Frequency [Hz]');ylabel('Phase [deg]'); xlim([7e-2 1e1]); %slide 6 figure subplot 211 loglog(Fr, abs(z12)./abs(z1));hold on; loglog(F_sl6,Amp_sl6,'*'); xlabel('Frequency [Hz]');ylabel('Amplitude [-]'); title('Z direction Top/bottom sensor'); xlim([7e-2 1e1]); subplot 212 semilogx(Fr, (180/pi).*angle(z12./abs(z1)));hold on; loglog(F_sl6_phase,Phase_sl6,'*'); xlabel('Frequency [Hz]');ylabel('Phase [deg]'); xlim([7e-2 1e1]);ylim([-180 180]); % %slide 6 % figure % loglog(Fr,abs(z2).^.5,Fr,abs(z1).^.5,Fr,abs(zsum).^.5,Fr,abs(zdelta).^.5) % xlabel('Frequency [Hz]');ylabel('Acceleration [m/s^2/rtHz]') % legend('Top sensor','Bottom sensor','Half sum','Half difference'); % title('Vertical') % xlim([7e-2 1e1]); %slide 9 figure subplot 211 loglog(Fr, abs(PHIwx)./abs(wx));hold on; loglog(F_sl9,Amp_sl9,'*'); xlabel('Frequency [Hz]');ylabel('Amplitude [-]'); title('X direction bottom/ground sensor'); xlim([7e-2 1e1]); ylim([0.01 10]); subplot 212 semilogx(Fr, (180/pi).*angle(PHIwx./abs(wx)));hold on; loglog(F_sl9,Phase_sl9,'*'); xlabel('Frequency [Hz]');ylabel('Phase [deg]'); xlim([7e-2 1e1]); % %slide 8 % figure % loglog(Fr,abs(wx).^.5,Fr,abs(x1).^.5,'-.',Fr,abs(x2).^.5,'.'); % grid on;xlabel('Frequency [Hz]'); % ylabel('ASD [m/s^2/rtHz]'); % xlim([7e-2 1e1]); % legend('Ground','Bottom sensor','Top sensor'); % title('Horizontal'); % % %slide 13 % figure % loglog(Fr,abs(wz).^.5,Fr,abs(z1).^.5,'-.',Fr,abs(z2).^.5,'.'); % grid on;xlabel('Frequency [Hz]'); % ylabel('ASD [m/s^2/rtHz]'); % xlim([7e-2 1e1]); % legend('Ground','Bottom sensor','Top sensor'); % title('Vertical'); %slide 14 figure subplot 211 loglog(Fr, abs(PHIwz)./abs(wz));hold on; loglog(F_sl14,Amp_sl14,'*'); xlabel('Frequency [Hz]');ylabel('Amplitude [-]'); title('Z direction bottom/ground sensor'); xlim([7e-2 1e1]); ylim([0.01 10]); subplot 212 semilogx(Fr, (180/pi).*angle(PHIwz./abs(wz)));hold on; loglog(F_sl14_phase,Phase_sl14,'*'); xlabel('Frequency [Hz]');ylabel('Phase [deg]'); xlim([7e-2 1e1]); %rotation figure loglog(Fr,abs(rot).^.5./((2*pi*Fr').^2),F_rot,Amp_rot,'*'); xlabel('Frequency [Hz]');ylabel('Rotation [rad/rtHz]') xlim([7e-2 1e1]); #+end_src ** Model Generation :noexport: #+begin_src matlab function [Fr, x1, z1, x2, z2, wx, wz, x12, z12, PHIwx, PHIwz,xsum,zsum,xdelta,zdelta,rot] = modelGeneration(m,I,k,h,ha,l,la,dampv,damph,kv,kh) %% generation of the state space model M = [m 0 0 0 m 0 0 0 I]; %Jacobian of the bottom sensor Js1 = [1 0 h/2 0 1 -l/2]; %Jacobian of the top sensor Js2 = [1 0 -h/2 0 1 0]; %Jacobian of the actuators Ja = [1 0 ha/2 %Left horizontal actuator %1 0 h/2 %Right horizontal actuator 0 1 -la/2 %Left vertical actuator 0 1 la/2]; %Right vertical actuator Jah = [1 0 ha/2]; Jav = [0 1 -la/2 %Left vertical actuator 0 1 la/2]; %Right vertical actuator Jta = Ja'; Jtah = Jah'; Jtav = Jav'; K = kv*Jtav*Jav + kh*Jtah*Jah; C = dampv*kv*Jtav*Jav+damph*kh*Jtah*Jah; E = [1 0 0 0 1 0 0 0 1]; %projecting ground motion in the directions of the legs AA = [zeros(3) eye(3) -M\K -M\C]; BB = [zeros(3,6) M\Jta M\(k*Jta*E)]; CC = [[Js1;Js2] zeros(4,3); zeros(2,6) (Js1+Js2)./2 zeros(2,3) (Js1-Js2)./2 zeros(2,3) (Js1-Js2)./(2*h) zeros(2,3)]; DD = [zeros(4,6) zeros(2,3) eye(2,3) zeros(6,6)]; system_dec = ss(AA,BB,CC,DD); %input = three actuators and three ground motions %output = the bottom sensor; the top sensor; the ground motion; the half %sum; the half difference; the rotation %% Injecting ground motion in the system to have the output Fr = logspace(-2,3,1e3); w=2*pi*Fr*1i; %fit of the ground motion data in m/s^2/rtHz Fr_ground_x = [0.07 0.1 0.15 0.3 0.7 0.8 0.9 1.2 5 10]; n_ground_x1 = [4e-7 4e-7 2e-6 1e-6 5e-7 5e-7 5e-7 1e-6 1e-5 3.5e-5]; Fr_ground_v = [0.07 0.08 0.1 0.11 0.12 0.15 0.25 0.6 0.8 1 1.2 1.6 2 6 10]; n_ground_v1 = [7e-7 7e-7 7e-7 1e-6 1.2e-6 1.5e-6 1e-6 9e-7 7e-7 7e-7 7e-7 1e-6 2e-6 1e-5 3e-5]; n_ground_x = interp1(Fr_ground_x,n_ground_x1,Fr,'linear'); n_ground_v = interp1(Fr_ground_v,n_ground_v1,Fr,'linear'); % figure % loglog(Fr,abs(n_ground_v),Fr_ground_v,n_ground_v1,'*'); % xlabel('Frequency [Hz]');ylabel('ASD [m/s^2 /rtHz]'); % return %converting into PSD n_ground_x = (n_ground_x).^2; n_ground_v = (n_ground_v).^2; %Injecting ground motion in the system and getting the outputs system_dec_f = (freqresp(system_dec,abs(w))); PHI = zeros(size(Fr,2),12,12); for p = 1:size(Fr,2) Sw=zeros(6,6); Iact = zeros(3,3); Sw(4,4) = n_ground_x(p); Sw(5,5) = n_ground_v(p); Sw(6,6) = n_ground_v(p); Sw(1:3,1:3) = Iact; PHI(p,:,:) = (system_dec_f(:,:,p))*Sw(:,:)*(system_dec_f(:,:,p))'; end x1 = PHI(:,1,1); z1 = PHI(:,2,2); x2 = PHI(:,3,3); z2 = PHI(:,4,4); wx = PHI(:,5,5); wz = PHI(:,6,6); x12 = PHI(:,1,3); z12 = PHI(:,2,4); PHIwx = PHI(:,1,5); PHIwz = PHI(:,2,6); xsum = PHI(:,7,7); zsum = PHI(:,8,8); xdelta = PHI(:,9,9); zdelta = PHI(:,10,10); rot = PHI(:,11,11); #+end_src * Simscape Model - Stewart Platform ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> #+end_src #+begin_src matlab :exports none :results silent :noweb yes <> #+end_src ** Jacobian First, the position of the "joints" (points of force application) are estimated and the Jacobian computed. #+begin_src matlab open('stewart_platform/drone_platform_jacobian.slx'); #+end_src #+begin_src matlab sim('drone_platform_jacobian'); #+end_src #+begin_src matlab Aa = [a1.Data(1,:); a2.Data(1,:); a3.Data(1,:); a4.Data(1,:); a5.Data(1,:); a6.Data(1,:)]'; Ab = [b1.Data(1,:); b2.Data(1,:); b3.Data(1,:); b4.Data(1,:); b5.Data(1,:); b6.Data(1,:)]'; As = (Ab - Aa)./vecnorm(Ab - Aa); l = vecnorm(Ab - Aa)'; J = [As' , cross(Ab, As)']; save('./jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J'); #+end_src ** Simulink #+begin_src matlab open('stewart_platform/drone_platform.slx'); #+end_src Definition of spring parameters #+begin_src matlab kx = 50; % [N/m] ky = 50; kz = 50; cx = 0.025; % [Nm/rad] cy = 0.025; cz = 0.025; #+end_src We load the Jacobian. #+begin_src matlab load('./jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J'); #+end_src The dynamics is identified from forces applied by each legs to the measured acceleration of the top platform. #+begin_src matlab %% Name of the Simulink File mdl = 'drone_platform'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/u'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Inertial Sensor'], 1, 'openoutput'); io_i = io_i + 1; G = linearize(mdl, io); G.InputName = {'F1', 'F2', 'F3', 'F4', 'F5', 'F6'}; G.OutputName = {'Ax', 'Ay', 'Az', 'Arx', 'Ary', 'Arz'}; #+end_src #+begin_src matlab :results output replace size(G) #+end_src #+RESULTS: : State-space model with 6 outputs, 6 inputs, and 12 states. Thanks to the Jacobian, we compute the transfer functions in the frame of the legs and in an inertial frame. #+begin_src matlab Gx = -G*inv(J'); Gx.InputName = {'Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'}; Gl = -J*G; Gl.OutputName = {'A1', 'A2', 'A3', 'A4', 'A5', 'A6'}; #+end_src #+begin_src matlab :exports none freqs = logspace(-1, 2, 1000); figure; ax1 = subplot(2, 1, 1); hold on; plot(freqs, abs(squeeze(freqresp(Gx(1, 1), freqs, 'Hz'))), 'DisplayName', '$A_x/F_x$'); plot(freqs, abs(squeeze(freqresp(Gx(2, 2), freqs, 'Hz'))), 'DisplayName', '$A_y/F_y$'); plot(freqs, abs(squeeze(freqresp(Gx(3, 3), freqs, 'Hz'))), 'DisplayName', '$A_z/F_z$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); legend('location', 'southeast'); ax2 = subplot(2, 1, 2); hold on; for i = 1:3 plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gx(i, i), freqs, 'Hz'))))); end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); ylim([-270, 90]); yticks([-360:90:360]); linkaxes([ax1,ax2],'x'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace exportFig('figs/stewart_platform_translations.pdf', 'width', 'full', 'height', 'full'); #+end_src #+name: fig:stewart_platform_translations #+caption: Stewart Platform Plant from forces applied by the legs to the acceleration of the platform #+RESULTS: [[file:figs/stewart_platform_translations.png]] #+begin_src matlab :exports none freqs = logspace(-1, 2, 1000); figure; ax1 = subplot(2, 1, 1); hold on; plot(freqs, abs(squeeze(freqresp(Gx(4, 4), freqs, 'Hz'))), 'DisplayName', '$A_{R_x}/M_x$'); plot(freqs, abs(squeeze(freqresp(Gx(5, 5), freqs, 'Hz'))), 'DisplayName', '$A_{R_y}/M_y$'); plot(freqs, abs(squeeze(freqresp(Gx(6, 6), freqs, 'Hz'))), 'DisplayName', '$A_{R_z}/M_z$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Amplitude [rad/(Nm)]'); set(gca, 'XTickLabel',[]); legend('location', 'southeast'); ax2 = subplot(2, 1, 2); hold on; for i = 4:6 plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gx(i, i), freqs, 'Hz'))))); end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); ylim([-270, 90]); yticks([-360:90:360]); linkaxes([ax1,ax2],'x'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace exportFig('figs/stewart_platform_rotations.pdf', 'width', 'full', 'height', 'full'); #+end_src #+name: fig:stewart_platform_rotations #+caption: Stewart Platform Plant from torques applied by the legs to the angular acceleration of the platform #+RESULTS: [[file:figs/stewart_platform_rotations.png]] #+begin_src matlab :exports none freqs = logspace(-1, 2, 1000); figure; ax1 = subplot(2, 1, 1); hold on; for i = 1:6 plot(freqs, abs(squeeze(freqresp(Gl(i, i), freqs, 'Hz')))); end for i = 1:5 for j = i+1:6 plot(freqs, abs(squeeze(freqresp(Gl(i, j), freqs, 'Hz'))), 'color', [0, 0, 0, 0.2]); end end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); ax2 = subplot(2, 1, 2); hold on; for i = 1:6 plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gl(i, i), freqs, 'Hz'))))); end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); ylim([-270, 90]); yticks([-360:90:360]); linkaxes([ax1,ax2],'x'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace exportFig('figs/stewart_platform_legs.pdf', 'width', 'full', 'height', 'full'); #+end_src #+name: fig:stewart_platform_legs #+caption: Stewart Platform Plant from forces applied by the legs to displacement of the legs #+RESULTS: [[file:figs/stewart_platform_legs.png]]