
%% Exported from Jupyter Notebook
% Run each section by placing your cursor in it and pressing Ctrl+Enter

%% Code Cell[1]:

set(0,'defaultTextInterpreter','latex');
clear; close all
save_figures = true;

%% Markdown Cell:
% Let's begin by defining the system parameters.

%% Code Cell[2]:

Kt_spec = 13.7; % oz-in/A ... torque constant from spec
Kv_spec = 10.2; % V/krpm ... voltage constant from spec
Tmax_spec = 2.82; % N-m ... max (stall) torque from spec
Omax_spec = 628; % rad/s ... max speed (no load) from spec 
N_oz = 0.278013851; % N/oz
m_in = 0.0254; % m/in
Kt_si = Kt_spec*N_oz*m_in; % N-m/A
rads_krpm = 1e3*2*pi/60; % (rad/s)/krpm
Kv_si = Kv_spec/rads_krpm; % V/(rad/s)
Jm = 56.5e-6; % kg-m^2 ... inertia of rotor
Bm = 16.9e-6; % N-m/s^2 ... motor damping coef
R = 1.6; % Ohm ... armature resistance
L = 4.1e-3; % H ... armature inductance
TF = Kv_si; % N-m/A ... trans ratio/motor constant

%% Markdown Cell:
% Let's investigate what happens in steady-state $\overline{\bm{x}}$.
% The system is stationary when $\dot{\bm{x}} = \bm{0}$ and $\bm{u} = \overline{\bm{u}}$ (stationary),\footnote{A stationary input $\overline{\bm{u}}$ is required for a stationary state if the input has any effect on the state; that is, if $B$ is nonzero.} so
% \begin{align}
%  \bm{0} &= A \overline{\bm{x}} + B \overline{\bm{u}} \Rightarrow \nonumber \\
%  \overline{\bm{x}} &= -A^{-1} B \overline{\bm{u}}.
% \end{align}
% Let's compute our steady-state solution for a constant voltage input $V_s(t) = \overline{V}$ and braking torque $T_s(t) = \overline{T}$.
% We use a symbolic approach to gain insight.

%% Code Cell[3]:

syms B_ J_ TF_ L_ R_ Vs_ Ts_ % using underscore for syms

a_ = [-B_/J_,TF_/J_;-TF_/L_,-R_/L_];
b_ = [0,-1/J_;1/L_,0];
u_ = [Vs_;Ts_];

M1_ = -inv(a_)*b_ % matrix -A^-1 B
den_ = TF_^2 + B_*R_; % common den
M2_ = M1_.*den_; % factor
xs_ = M1_*u_ % full ss sol
xs_2_ = M2_*u_; % naughty factorless ss sol

%% Code Cell[4]:

eig(a_)

%% Markdown Cell:
% A little more human-readably, using the fact that $\Omega_2 = \Omega_J$ and $i_1 = i_L$, and using bars to denote steady-state values,
% \begin{align}
% \overline{\Omega_2} &= \dfrac{1}{TF^2 + B_m R}(TF\, \overline{V_s} - R \overline{T_s}) \\
% \overline{i_1} &= \dfrac{1}{TF^2 + B_m R}(B_m \overline{V_s} + TF\, \overline{T_s})
% \end{align}

%% Markdown Cell:
% Let's focus on the first of these, the relationship between $\overline{\Omega_2}$ and $\overline{T_s}$.
% For given $\overline{V_s}$, there is a linearly decreasing relationship between $\overline{\Omega_2}$ and $\overline{T_s}$. 
% This is precisely the *motor curve*.
% But it's one of a few curves plotted versus $\overline{T_s}$.
% Other common curves are current $\overline{i_1}$, mechanical braking power $\mathcal{P}_\text{brk} = \overline{T_s} \overline{\Omega_s}$, and efficiency $\varepsilon$.
% The efficiency is defined as the ratio of the braking power to the voltage source power $\mathcal{P}_\text{src} = \overline{I_s} \overline{V_s}$; i.e. 
% \begin{align}
%     \varepsilon = \mathcal{P}_\text{brk}/\mathcal{P}_\text{src}.
% \end{align}
% 
% We already have expressions for $\overline{\Omega_2}$ and $\overline{i_1}$ in terms of $\overline{T_s}$, but we must still derive such expressions for $\mathcal{P}_\text{brk}$ and $\varepsilon$.
% For $\mathcal{P}_\text{brk}$, we must express $\overline{\Omega_s}$ in terms of known quantities.
% From the linear graph, it is obvious that $\overline{\Omega_s} = \overline{\Omega_2}$.
% Therefore,
% \begin{align}
%     \mathcal{P}_\text{brk} = \overline{T_s} \overline{\Omega_2}.
% \end{align}
% Now for $\varepsilon$.
% We have the unknown source current $\overline{I_s}$.
% However, from the linear graph, it is obvious that $\overline{I_s} = \overline{i_1}$.
% Therefore,
% \begin{align}
%     \varepsilon = \frac{\overline{T_s}\overline{\Omega_2}}{\overline{i_1}\overline{V_s}}.
% \end{align}
% 
% Let's compute these quantities for our parameters.

%% Code Cell[5]:

Vs = 60; % V ... max used, which is common
Tmax = TF/R*Vs; % N-m ... occurs when Omega_J = 0
Ts_a = linspace(0,Tmax,180); % N-m ... braking torques
O2_a = 1/(TF^2 + Bm*R)*(TF*Vs-R*Ts_a); % rad/s ... ss speed
i1_a = 1/(TF^2 + Bm*R)*(Bm*Vs+TF*Ts_a);
Pbrk_a = Ts_a.*O2_a; % W ... braking power
eff_a = Pbrk_a./(i1_a*Vs);

%% Markdown Cell:
% Now let's plot them!
% The output is shown in \autoref{fig:embedded_motor_model_motor_curves}.

%% Code Cell[6]:

data(1).X = Ts_a;
data(1).Y = O2_a; 
data(1).YLabel = '$\overline{\Omega_2}$ (rad/s)';
data(2).Y = i1_a; 
data(2).YLabel = '$\overline{i_1}$ (A)';
data(3).Y = Pbrk_a; 
data(3).YLabel = '$\mathcal{P}_{\mathrm{brk}} (W)$';
data(4).Y = eff_a; 
data(4).YLabel = '$\varepsilon$';
data(1).XLabel = '$\overline{T_s}$';
data(1).LineWidth = 1.5;
data(2).LineWidth = 1.5;
data(3).LineWidth = 1.5;
data(4).LineWidth = 1.5;
h = figure;
[hax, hlines, data] = plotyn(data);
hgsave(h,'figures/temp');

%% Code Cell[7]:

if save_figures
    fn = 'embedded_motor_model_motor_curves.tex';
    h = hgload('figures/temp');
    cleanfigure;
    matlab2tikz(['figures/',fn],...
        'figurehandle',h,...
        'width','.9\linewidth',...
        'height','.6\linewidth',...
        'extraAxisOptions',...
        [...
            'every axis plot/.append style={thick},',...
            'scaled ticks=false,'
        ],...
        'showInfo', false...
    );
end

%% Markdown Cell:
% There are some key quantities that can be read from the graph and found analytically.
% The most important are the maximum speed $\overline{\Omega_2}_\mathrm{max}$, which occurs at zero torque, and maximum torque $\overline{T_s}_\mathrm{max}$, which occurs at zero speed.
% Another is that the maximum mechanical power (output) occurs at $\overline{T_s}_\mathrm{max}/2$.
% Finally, the maximum effieciency occurs at relatively low torque and high speed, which is typical for the following reason: the two energy-dissipative elements, the resistor and the damper, trade-off as being the dominant effect at the peak, and the resistor tends to dominate.
% That is, at high speed/voltage and low torque/current, the damper dominates dissipation;
% at low speed/voltage and high torque/current, the resistor dominates dissipation.
% It is very common for a motor's resistance to dominate the damping, as in our case.
% 
% Let's examine the maximum speed and torque.

%% Code Cell[8]:

Omax = O2_a(1) % rad/s ... occurs when T_s = 0
Tmax % N-m ... already computed and occurs when Omega_2 = 0

%% Markdown Cell:
% Comparing these to the values given in the spec sheet, we see we're pretty good, but there's a bit of a discrepency in the max torque.

%% Code Cell[10]:

Omax_spec
Tmax_spec
disp(sprintf('percent error for speed: %0.3g',...
    (Omax-Omax_spec)/Omax_spec*100))
disp(sprintf('percent error for torque: %0.3g',...
    (Tmax-Tmax_spec)/Tmax_spec*100))

%% Markdown Cell:
% We should investigate further, but what we will find is that these values are fairly sensitive to $TF$, $B$, and $R$.
% In our case, it is likely that the given value for $R$ is a bit low.
% It is given as $1.6\ \Omega$, but it is probably closer to $2\ \Omega$.
% However, the datasheet for this motor was not clear about whether the maximum speed and torque values were derived from a full motor curve fit or if they were the only points measured.
% The former is best for estimating dynamic model parameters like $R$ and $B$, but the latter is occasionally sufficient.

%% Code Cell[11]:

%TODO: plot all the steady-state outputs, perhaps in a stack or all on the same plot
