diff --git a/Figures/G_ws_pz.pdf b/Figures/G_ws_pz.pdf new file mode 100644 index 0000000..f4aa94d Binary files /dev/null and b/Figures/G_ws_pz.pdf differ diff --git a/Figures/G_ws_pz.png b/Figures/G_ws_pz.png new file mode 100644 index 0000000..aba5764 Binary files /dev/null and b/Figures/G_ws_pz.png differ diff --git a/Figures/G_ws_pz.svg b/Figures/G_ws_pz.svg new file mode 100644 index 0000000..0d21b18 Binary files /dev/null and b/Figures/G_ws_pz.svg differ diff --git a/Figures/G_ws_vc.pdf b/Figures/G_ws_vc.pdf new file mode 100644 index 0000000..eda05d0 Binary files /dev/null and b/Figures/G_ws_vc.pdf differ diff --git a/Figures/G_ws_vc.png b/Figures/G_ws_vc.png new file mode 100644 index 0000000..f730aa8 Binary files /dev/null and b/Figures/G_ws_vc.png differ diff --git a/Figures/G_ws_vc.svg b/Figures/G_ws_vc.svg new file mode 100644 index 0000000..80f922c Binary files /dev/null and b/Figures/G_ws_vc.svg differ diff --git a/Figures/Gc_ws_pz.pdf b/Figures/Gc_ws_pz.pdf new file mode 100644 index 0000000..cb97deb Binary files /dev/null and b/Figures/Gc_ws_pz.pdf differ diff --git a/Figures/Gc_ws_pz.png b/Figures/Gc_ws_pz.png new file mode 100644 index 0000000..60f3dae Binary files /dev/null and b/Figures/Gc_ws_pz.png differ diff --git a/Figures/Gc_ws_pz.svg b/Figures/Gc_ws_pz.svg new file mode 100644 index 0000000..8245d76 Binary files /dev/null and b/Figures/Gc_ws_pz.svg differ diff --git a/Figures/Gc_ws_vc.pdf b/Figures/Gc_ws_vc.pdf new file mode 100644 index 0000000..4376991 Binary files /dev/null and b/Figures/Gc_ws_vc.pdf differ diff --git a/Figures/Gc_ws_vc.png b/Figures/Gc_ws_vc.png new file mode 100644 index 0000000..eac4120 Binary files /dev/null and b/Figures/Gc_ws_vc.png differ diff --git a/Figures/Gc_ws_vc.svg b/Figures/Gc_ws_vc.svg new file mode 100644 index 0000000..1ae0a2c Binary files /dev/null and b/Figures/Gc_ws_vc.svg differ diff --git a/Figures/poles_w_pz.pdf b/Figures/poles_w_pz.pdf new file mode 100644 index 0000000..2f4e221 Binary files /dev/null and b/Figures/poles_w_pz.pdf differ diff --git a/Figures/poles_w_pz.png b/Figures/poles_w_pz.png new file mode 100644 index 0000000..db1f7eb Binary files /dev/null and b/Figures/poles_w_pz.png differ diff --git a/Figures/poles_w_pz.svg b/Figures/poles_w_pz.svg new file mode 100644 index 0000000..61b70dd Binary files /dev/null and b/Figures/poles_w_pz.svg differ diff --git a/Figures/poles_w_vc.pdf b/Figures/poles_w_vc.pdf new file mode 100644 index 0000000..daceb2c Binary files /dev/null and b/Figures/poles_w_vc.pdf differ diff --git a/Figures/poles_w_vc.png b/Figures/poles_w_vc.png new file mode 100644 index 0000000..b11dd02 Binary files /dev/null and b/Figures/poles_w_vc.png differ diff --git a/Figures/poles_w_vc.svg b/Figures/poles_w_vc.svg new file mode 100644 index 0000000..ddf24fb Binary files /dev/null and b/Figures/poles_w_vc.svg differ diff --git a/mat/parameters.mat b/mat/parameters.mat index 41134f2..1c28b4f 100644 Binary files a/mat/parameters.mat and b/mat/parameters.mat differ diff --git a/rotating_frame.html b/rotating_frame.html index a2f1cd0..aa7c36c 100644 Binary files a/rotating_frame.html and b/rotating_frame.html differ diff --git a/rotating_frame.org b/rotating_frame.org index b1b8ab9..9eb1755 100644 --- a/rotating_frame.org +++ b/rotating_frame.org @@ -33,7 +33,7 @@ In section [[sec:simscape]], the analytical study will be validated using a mult Finally, in section [[sec:control]], the control strategies are implemented using Simulink and Simscape and compared. -* System +* System Description and Analysis :PROPERTIES: :HEADER-ARGS:matlab+: :tangle system_numerical_analysis.m :END: @@ -290,6 +290,7 @@ Similarly we can obtain $d_v$ function of $F_u$ and $F_v$: The two previous equations can be written in a matrix form: #+begin_important +#+NAME: eq:coupledplant \begin{equation} \begin{bmatrix} d_u \\ d_v \end{bmatrix} = \frac{1}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} @@ -308,8 +309,7 @@ We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \ome - with the light sample (figure [[fig:coupling_light]]). - with the heavy sample (figure [[fig:coupling_heavy]]). -#+HEADER: :exports none :results silent -#+begin_src matlab +#+begin_src matlab :exports none :results silent f = logspace(-1, 3, 1000); figure; @@ -335,8 +335,7 @@ We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \ome #+RESULTS: [[file:Figures/coupling_light.png]] -#+HEADER: :exports none :results silent -#+begin_src matlab +#+begin_src matlab :exports none :results silent figure; hold on; plot(f, abs(-mheavy*(2*pi*f).^2 + kvc - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Voice Coil') @@ -405,6 +404,342 @@ From this analysis, we can determine the lowest practical stiffness that is poss |-------------+-------+-------| | k min [N/m] | 2199 | 89 | +** Effect of rotation speed on the plant +As shown in equation [[eq:coupledplant]], the plant changes with the rotation speed $\omega_0$. + +Then, we compute the bode plot of the direct term and coupling term for multiple rotating speed. + +Then we compare the result between voice coil and piezoelectric actuators. + +*** Voice coil actuator + +#+begin_src matlab :exports none :results silent + m = mlight; + k = kvc; +#+end_src + +#+begin_src matlab :exports none :results silent + ws = [0, 1, 10, 60]*2*pi/60; % [rmp] + + Gs = {zeros(1, length(ws))}; + Gcs = {zeros(1, length(ws))}; + + for i = 1:length(ws) + w = ws(i); + Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)}; + Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)}; + end +#+end_src + +#+begin_src matlab :exports none :results silent + freqs = logspace(-2, 1, 1000); + + figure; + ax1 = subaxis(2,1,1); + hold on; + for i = 1:length(ws) + plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi)); + end + hold off; + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [m/N]'); + legend('Location', 'southwest'); + + ax2 = subaxis(2,1,2); + hold on; + for i = 1:length(ws) + plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz')))); + end + hold off; + yticks(-180:90:180); + ylim([-180 180]); + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + linkaxes([ax1,ax2],'x'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/G_ws_vc.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:G_ws_vc +#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Voice coil +#+RESULTS: +[[file:Figures/G_ws_vc.png]] + +#+begin_src matlab :exports none :results silent + freqs = logspace(-2, 1, 1000); + + figure; + ax1 = subaxis(2,1,1); + hold on; + for i = 1:length(ws) + plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi)); + end + hold off; + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [m/N]'); + legend('Location', 'southwest'); + + ax2 = subaxis(2,1,2); + hold on; + for i = 1:length(ws) + plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz')))); + end + hold off; + yticks(-180:90:180); + ylim([-180 180]); + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + linkaxes([ax1,ax2],'x'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/Gc_ws_vc.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:Gc_ws_vc +#+CAPTION: Bode plot of the coupling transfer function term (from $F_u$ to $D_v$) for multiple rotation speed - Voice coil +#+RESULTS: +[[file:Figures/Gc_ws_vc.png]] + +*** Piezoelectric actuator + +#+begin_src matlab :exports none :results silent + m = mlight; + k = kpz; +#+end_src + +#+begin_src matlab :exports none :results silent + ws = [0, 1, 10, 60]*2*pi/60; % [rmp] + + Gs = {zeros(1, length(ws))}; + Gcs = {zeros(1, length(ws))}; + + for i = 1:length(ws) + w = ws(i); + Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)}; + Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)}; + end +#+end_src + +#+begin_src matlab :exports none :results silent + freqs = logspace(2, 3, 1000); + + figure; + ax1 = subaxis(2,1,1); + hold on; + for i = 1:length(ws) + plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi)); + end + hold off; + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [m/N]'); + legend('Location', 'southwest'); + + ax2 = subaxis(2,1,2); + hold on; + for i = 1:length(ws) + plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz')))); + end + hold off; + yticks(-180:90:180); + ylim([-180 180]); + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + linkaxes([ax1,ax2],'x'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/G_ws_pz.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:G_ws_pz +#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Piezoelectric actuator +#+RESULTS: +[[file:Figures/G_ws_pz.png]] + +#+begin_src matlab :exports none :results silent + figure; + ax1 = subaxis(2,1,1); + hold on; + for i = 1:length(ws) + plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi)); + end + hold off; + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [m/N]'); + legend('Location', 'southwest'); + + ax2 = subaxis(2,1,2); + hold on; + for i = 1:length(ws) + plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz')))); + end + hold off; + yticks(-180:90:180); + ylim([-180 180]); + xlim([freqs(1), freqs(end)]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + linkaxes([ax1,ax2],'x'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/Gc_ws_pz.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:Gc_ws_pz +#+CAPTION: Bode plot of the coupling transfer function term (from $F_u$ to $D_v$) for multiple rotation speed - Piezoelectric actuator +#+RESULTS: +[[file:Figures/Gc_ws_pz.png]] + +*** Analysis +When the rotation speed is null, the coupling terms are equal to zero and the diagonal terms corresponds to one degree of freedom mass spring system. + +When the rotation speed in not null, the resonance frequency is duplicated into two pairs of complex conjugate poles. + +As the rotation speed increases, one of the two resonant frequency goes to lower frequencies as the other one goes to higher frequencies. + +The poles of the coupling terms are the same as the poles of the diagonal terms. The magnitude of the coupling terms are increasing with the rotation speed. + +As shown in the previous figures, the system with voice coil is much more sensitive to rotation speed. + +*** Campbell diagram + +The poles of the system are computed for multiple values of the rotation frequency. + +#+begin_src matlab :results silent :exports code + m = mlight; + k = kvc; + c = 0.1*sqrt(k*m); + + ws = linspace(0, 10, 100); % [rad/s] + + polesvc = zeros(2, length(ws)); + + for i = 1:length(ws) + polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2)); + polesvc(:, i) = sort(polei(imag(polei) > 0)); + end +#+end_src + +#+begin_src matlab :results silent :exports code + m = mlight; + k = kpz; + c = 0.1*sqrt(k*m); + + ws = linspace(0, 1000, 100); % [rad/s] + + polespz = zeros(2, length(ws)); + + for i = 1:length(ws) + polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2)); + polespz(:, i) = sort(polei(imag(polei) > 0)); + end +#+end_src + +We then plot the real and imaginary part of the poles as a function of the rotation frequency (figures [[fig:poles_w_vc]] and [[fig:poles_w_pz]]). + +When the real part of one pole becomes positive, the system goes unstable. + +For the voice coil (figure [[fig:poles_w_vc]]), the system is unstable when the rotation speed is above 5 rad/s. The real and imaginary part of the poles of the system with piezoelectric actuators are changing much less (figure [[fig:poles_w_pz]]). + +#+begin_src matlab :results silent :exports none + figure; + % Amplitude + ax1 = subplot(1,2,1); + hold on; + plot(ws, real(polesvc(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, real(polesvc(2, :))) + plot(ws, zeros(size(ws)), 'k--') + hold off; + xlabel('Rotation Frequency [rad/s]'); + ylabel('Pole Real Part'); + ax2 = subplot(1,2,2); + hold on; + plot(ws, imag(polesvc(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, -imag(polesvc(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, imag(polesvc(2, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, -imag(polesvc(2, :))) + hold off; + xlabel('Rotation Frequency [rad/s]'); + ylabel('Pole Imaginary Part'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/poles_w_vc.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:poles_w_vc +#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample +#+RESULTS: +[[file:Figures/poles_w_vc.png]] + + +#+begin_src matlab :results silent :exports none + figure; + % Amplitude + ax1 = subplot(1,2,1); + hold on; + plot(ws, real(polespz(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, real(polespz(2, :))) + plot(ws, zeros(size(ws)), 'k--') + hold off; + xlabel('Rotation Frequency [rad/s]'); + ylabel('Pole Real Part'); + ax2 = subplot(1,2,2); + hold on; + plot(ws, imag(polespz(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, -imag(polespz(1, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, imag(polespz(2, :))) + set(gca,'ColorOrderIndex',1) + plot(ws, -imag(polespz(2, :))) + hold off; + xlabel('Rotation Frequency [rad/s]'); + ylabel('Pole Imaginary Part'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/poles_w_pz.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + +#+NAME: fig:poles_w_pz +#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample +#+RESULTS: +[[file:Figures/poles_w_pz.png]] + + * Control Strategies <> ** Measurement in the fixed reference frame @@ -445,6 +780,7 @@ The loop gain is $L = G K$. :END: <> +** Initialize #+begin_src matlab :exports none :results silent :noweb yes <> load('./mat/parameters.mat'); @@ -455,6 +791,7 @@ The loop gain is $L = G K$. #+end_src ** Parameter for the Simscape simulations +First we define the parameters that must be defined in order to run the Simscape simulation. #+begin_src matlab :exports code :results silent w = 2*pi; % Rotation speed [rad/s] @@ -467,6 +804,7 @@ The loop gain is $L = G K$. cTuv = 0; % Damping of the Tuv stage [N/(m/s)] #+end_src +Then, we defined parameters that will be used in the following analysis. #+begin_src matlab :exports code :results silent mlight = 5; % Mass for light sample [kg] mheavy = 55; % Mass for heavy sample [kg] @@ -478,11 +816,15 @@ The loop gain is $L = G K$. kpz = 1e8; % Piezo Stiffness [N/m] d = 0.01; % Maximum excentricity from rotational axis [m] + + freqs = logspace(-2, 3, 1000); % Frequency vector for analysis [Hz] #+end_src ** Identification in the rotating referenced frame +We initialize the inputs and outputs of the system to identify: +- Inputs: $f_u$ and $f_v$ +- Outputs: $d_u$ and $d_v$ -We initialize the inputs and outputs of the system to identify. #+begin_src matlab :exports code :results silent %% Options for Linearized options = linearizeOptions; @@ -531,11 +873,11 @@ Then we identify the system with an heavy mass and low speed. Gvc_heavy.OutputName = {'Du', 'Dv'}; #+end_src -Finally, we plot the coupling ratio in both case (figure [[fig:coupling_ration_light_heavy]]). +** Coupling ratio between $f_{uv}$ and $d_{uv}$ + +From the previous identification, we plot the coupling ratio in both case (figure [[fig:coupling_ration_light_heavy]]). We obtain the same result than the analytical case (figures [[fig:coupling_light]] and [[fig:coupling_heavy]]). #+begin_src matlab :results silent :exports none - freqs = logspace(-2, 3, 1000); - figure; hold on; plot(freqs, abs(squeeze(freqresp(Gvc_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_light('Dv', 'Fu'), freqs, 'Hz')))); @@ -560,6 +902,108 @@ We obtain the same result than the analytical case (figures [[fig:coupling_light #+RESULTS: [[file:Figures/coupling_ration_light_heavy.png]] +** Plant Control + +The goal is the study control problems due to the coupling that appears because of the rotation. + +#+begin_src matlab :exports none :results silent + %% Options for Linearized + options = linearizeOptions; + options.SampleTime = 0; + + %% Name of the Simulink File + mdl = 'rotating_frame'; + + %% Input/Output definition + io(1) = linio([mdl, '/fu'], 1, 'input'); + io(2) = linio([mdl, '/fv'], 1, 'input'); + + io(3) = linio([mdl, '/du'], 1, 'output'); + io(4) = linio([mdl, '/dv'], 1, 'output'); +#+end_src + +First, we identify the system when the rotation speed is null and then when the rotation speed is equal to 60rpm. + +The actuators are voice coil with some damping. + +#+begin_src matlab :exports none :results silent + w = 0; % Rotation speed [rad/s] + m = mlight; % mass of the sample [kg] + kTuv = kvc; + % cTuv = 0.1*sqrt(kTuv*m); + cTuv = 0; + + G = linearize(mdl, io, 0.1); + G.InputName = {'Fu', 'Fv'}; + G.OutputName = {'Du', 'Dv'}; +#+end_src + +#+begin_src matlab :exports none :results silent + w = 0.1; % Rotation speed [rad/s] + m = mlight; % mass of the sample [kg] + kTuv = kvc; + % cTuv = 0.1*sqrt(kTuv*m); + cTuv = 0; + + Gt = linearize(mdl, io, 0.1); + Gt.InputName = {'Fu', 'Fv'}; + Gt.OutputName = {'Du', 'Dv'}; +#+end_src + +#+begin_src matlab :exports none :results silent + figure; + bode(G, 'r-', Gt, 'b--') + legend({'G - w = 0', 'G - w = 60rpm'}, 'Location', 'southwest'); +#+end_src + +#+HEADER: :tangle no :exports results :results file :noweb yes +#+HEADER: :var filepath="Figures/coupling_ration_light_heavy.png" :var figsize="wide-tall" +#+begin_src matlab + <> +#+end_src + + +#+begin_src matlab :exports none :results silent + figure; + % Amplitude + ax1 = subaxis(2,1,1); + hold on; + plot(freqs, abs(squeeze(freqresp(Gvc_light('du', 'fu'), freqs, 'Hz'))), '-'); + set(gca,'xscale','log'); set(gca,'yscale','log'); + ylabel('Amplitude [m/N]'); + set(gca, 'XTickLabel',[]); + hold off; + % Phase + ax2 = subaxis(2,1,2); + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gvc_light('du', 'fu'), freqs, 'Hz')))), '-'); + set(gca,'xscale','log'); + yticks(-180:180:180); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); +#+end_src + +#+begin_src matlab :exports none :results silent + figure; + % Amplitude + ax1 = subaxis(2,1,1); + hold on; + plot(freqs, abs(squeeze(freqresp(-Gvc_light('dv', 'fv'), freqs, 'Hz'))), '-'); + set(gca,'xscale','log'); set(gca,'yscale','log'); + ylabel('Amplitude [m/N]'); + set(gca, 'XTickLabel',[]); + hold off; + % Phase + ax2 = subaxis(2,1,2); + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(-Gvc_light('dv', 'fv'), freqs, 'Hz')))), '-'); + set(gca,'xscale','log'); + yticks(-180:180:180); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); +#+end_src