\documentclass[a4paper,14pt]{extarticle}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[russian]{babel}

\usepackage[a4paper,margin=2.5cm]{geometry}
\usepackage{array}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage{amsmath,amssymb}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{float}
\usepackage{titlesec}
\usepackage{indentfirst}
\usepackage{newtxtext,newtxmath}

\linespread{1.3}
\frenchspacing
\setlength{\parindent}{1.25cm}
\setlength{\parskip}{0pt}

\usepackage[T2A]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[russian]{babel}

\usepackage{xcolor}
\usepackage{matlab-prettifier}

% Общие настройки листингов
\lstset{
  style=Matlab-editor,      % стиль из matlab-prettifier
  basicstyle=\mlttfamily\small, % моноширинный шрифт поменьше
  numbers=left,
  numberstyle=\tiny,
  stepnumber=1,
  numbersep=5pt,
  frame=single,
  breaklines=true,
  breakatwhitespace=false,
  captionpos=b,
  keepspaces=true,
  columns=fullflexible,
}

\lstset{style=matlabstyle}

\begin{document}

\tableofcontents
\newpage

% ---------- Введение ----------
\section*{Введение}
\addcontentsline{toc}{section}{Введение}

Целью лабораторной работы является исследование метода модального управления
дискретными системами автоматического управления и влияния размещения полюсов
замкнутой системы на качество переходного процесса.

В ходе работы требуется:
\begin{itemize}
  \item задать модель объекта управления в пространстве состояний;
  \item выполнить её дискретизацию для нескольких значений периода дискретизации $T_0$;
  \item проверить управляемость дискретной системы;
  \item синтезировать регулятор по состоянию методом модального управления
        для двух вариантов размещения полюсов: с перерегулированием порядка $30\%$
        и с апериодическим переходным процессом;
  \item подобрать масштабирующий коэффициент $k_0$,
        обеспечивающий нулевую установившуюся ошибку на ступенчатом воздействии;
  \item реализовать структуру замкнутой системы в среде Simulink без применения
        блока \texttt{State-Space} и исследовать переходные процепссы.
\end{itemize}

% ---------- 1. Исходные данные ----------
\section{Исходные данные}

Объект управления задан в виде
\begin{equation}
\dot x(t) = A x(t) + B u(t), \qquad y(t) = C x(t),
\end{equation}
где $x(t)\in\mathbb{R}^3$~--- вектор состояния,
$u(t)$~--- управляющее воздействие, $y(t)$~--- выход системы.

Для варианта~1 матрицы имеют вид
\begin{equation}
A=
\begin{bmatrix}
-39   & -71   & -23 \\
 23   &  40   &  13 \\
-6{,}6 & -10{,}7 & -3{,}3
\end{bmatrix},\quad
B=
\begin{bmatrix}
 0{,}228 \\
-0{,}127 \\
 0{,}038
\end{bmatrix},\quad
C=
\begin{bmatrix}
28000 & 84000 & 112000
\end{bmatrix}.
\end{equation}

Периоды дискретизации, используемые в работе:
\[
T_0 \in \{0{,}1;\ 1;\ 10\}\ \text{с}.
\]

% ---------- 2. Дискретизация и управляемость ----------
\section{Дискретизация модели и проверка управляемости}

При дискретизации по нулевому порядку удержания ($ZOH$) непрерывная модель
переходит в дискретную:
\begin{equation}
x_{k+1} = A_d x_k + B_d u_k,\qquad y_k = C_d x_k,
\end{equation}
где
\begin{equation}
A_d = e^{A T_0},\qquad
B_d = \int_{0}^{T_0} e^{A\tau} B\, d\tau,\qquad
C_d = C.
\end{equation}

Расчёты выполнены в среде \texttt{MATLAB}. Получены следующие матрицы.

\subsection*{Период дискретизации $T_0=0{,}1$~с}
\[
A_d^{(0{,}1)} =
\begin{bmatrix}
 -2{,}3833 & -5{,}7090 & -1{,}8796 \\
  1{,}8225 &  3{,}8935 &  0{,}9575 \\
 -0{,}4563 & -0{,}6412 &  0{,}8052
\end{bmatrix},\quad
B_d^{(0{,}1)} =
\begin{bmatrix}
 0{,}017508 \\
-0{,}0087103 \\
 0{,}0023217
\end{bmatrix}.
\]

\subsection*{Период дискретизации $T_0=1$~с}
\[
A_d^{(1)} =
\begin{bmatrix}
 -1{,}2509 & -2{,}5863 & -0{,}9697 \\
  0{,}7203 &  1{,}2968 &  0{,}0901 \\
  0{,}0844 &  0{,}5675 &  1{,}4380
\end{bmatrix},\quad
B_d^{(1)} =
\begin{bmatrix}
 0{,}0048087 \\
-0{,}00085652 \\
 0{,}0041248
\end{bmatrix}.
\]

\subsection*{Период дискретизации $T_0=10$~с}
\[
A_d^{(10)} =
\begin{bmatrix}
 -0{,}39861 & -1{,}1924 & -1{,}7848 \\
 -0{,}85955 & -2{,}5711 & -3{,}8485 \\
  3{,}3336  &  9{,}9717 & 14{,}9260
\end{bmatrix},\quad
B_d^{(10)} =
\begin{bmatrix}
-0{,}029755 \\
-0{,}052071 \\
 0{,}22143
\end{bmatrix}.
\]

\subsection*{Проверка управляемости}

Матрица управляемости дискретной системы:
\[
W_c = \bigl[B_d,\ A_d B_d,\ A_d^2 B_d\bigr].
\]

Для всех трёх значений периода дискретизации ранги матриц управляемости равны
порядку системы:
\[
\rank W_c^{(0{,}1)} = \rank W_c^{(1)} = \rank W_c^{(10)} = 3.
\]
Следовательно, дискретные модели объекта полностью управляемы, что позволяет
использовать метод модального управления.s

% ---------- 3. Синтез модального регулятора ----------
\section{Синтез модального регулятора по состоянию}

Примем закон управления по состоянию
\begin{equation}
u_k = -K x_k + k_0 g_k,
\end{equation}
где $K = \begin{bmatrix}k_1 & k_2 & k_3\end{bmatrix}$~--- вектор обратной
связи по состоянию, $k_0$~--- масштабирующий коэффициент, $g_k$~---
задающее воздействие (ступенька единичной амплитуды).

Матрица замкнутой системы:
\begin{equation}
A_{\text{cl}} = A_d - B_d K.
\end{equation}

Собственные значения матрицы $A_{\text{cl}}$ задаются заранее согласно
требуемому характеру переходного процесса.

\subsection{Выбор желаемых полюсов}

Для варианта с перерегулированием $\approx30\%$ используется связь
между перерегулированием $M_p$ и коэффициентом затухания $\zeta$:
\begin{equation}
M_p = e^{-\pi \zeta / \sqrt{1-\zeta^2}} \approx 0{,}3,
\end{equation}
откуда $\zeta \approx 0{,}36$.

В непрерывной области выбирается доминирующая комплексно-сопряжённая пара
полюсов
\[
s_{1,2} = -\sigma \pm j \omega_d,
\]
а также дополнительный действительный полюс $s_3=-5$.
Для малых и средних значений периода дискретизации ($T_0=0{,}1$ и $T_0=1$)
желаемые дискретные полюса получены экспоненциальным отображением
\[
z_i = e^{s_i T_0}.
\]

Для крупного шага $T_0=10$ устойчивые дискретные полюса заданы непосредственно
в $z$-области:
\[
z_{30} = \{0{,}7\pm0{,}3j;\ 0{,}4\},
\qquad
z_0 = \{0{,}7;\ 0{,}5;\ 0{,}3\}.
\]

\subsection{Расчёт коэффициентов $K$}

Вектор обратной связи $K$ находится методом Аккерманна
(\texttt{place} в \texttt{MATLAB}) для каждого периода дискретизации и
каждого варианта размещения полюсов.

Результаты сведены в таблицы
\ref{tab:K30} и~\ref{tab:K0}.

\begin{table}[H]
\centering
\caption{Коэффициенты $K$ для варианта с перерегулированием $\approx30\%$}
\label{tab:K30}
\begin{tabular}{cccc}
\toprule
$T_0$, c & $k_1$ & $k_2$ & $k_3$ \\
\midrule
0{,}1 &  209{,}9   &  592{,}6   &  868{,}0   \\
1     &   72{,}65  &  208{,}0   &  289{,}4   \\
10    & 598900{,}0 & 390900{,}0 & 172500{,}0 \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[H]
\centering
\caption{Коэффициенты $K$ для апериодического варианта}
\label{tab:K0}
\begin{tabular}{cccc}
\toprule
$T_0$, c & $k_1$ & $k_2$ & $k_3$ \\
\midrule
0{,}1 &  -54{,}92  & -108{,}0  &  45{,}88  \\
1     &   69{,}46  &  195{,}9  & 270{,}1  \\
10    & 271100{,}0 & 176900{,}0 & 78080{,}0 \\
\bottomrule
\end{tabular}
\end{table}

Собственные значения матриц $A_{\text{cl}}$ во всех случаях совпадают с
заданными полюсами, что подтверждает корректность синтеза регулятора.

% ---------- 4. Масштабирующий коэффициент k0 ----------
\section{Определение масштабирующего коэффициента $k_0$}

Для замкнутой системы со ступенчатым входом $g_k=1$ установившееся значение
выхода определяется выражением
\begin{equation}
y_\infty = C_d (I - A_{\text{cl}})^{-1} B_d\, k_0,
\qquad
A_{\text{cl}} = A_d - B_d K.
\end{equation}

Требуя $y_\infty = 1$, получаем формулу для масштабирующего коэффициента:
\begin{equation}
k_0 = \frac{1}{C_d (I - A_{\text{cl}})^{-1} B_d}.
\end{equation}

Численные значения $k_0$ для всех режимов приведены в таблице~\ref{tab:k0}.

\begin{table}[H]
\centering
\caption{Масштабирующий коэффициент $k_0$}
\label{tab:k0}
\begin{tabular}{ccc}
\toprule
$T_0$, c & $k_0$ (30\% перерегулир.) & $k_0$ (апериодич.) \\
\midrule
0{,}1 & $7{,}547\cdot 10^{-3}$ & $7{,}313\cdot 10^{-4}$ \\
1     & $2{,}066\cdot 10^{-3}$ & $1{,}88\cdot 10^{-3}$  \\
10    & $5{,}669\cdot 10^{-6}$ & $5{,}511\cdot 10^{-6}$ \\
\bottomrule
\end{tabular}
\end{table}

Подбор $k_0$ обеспечивает нулевую установившуюся ошибку по ступенчатому
задающему воздействию во всех рассмотренных случаях.s

% ---------- 5. Моделирование и графики ----------
\section{Моделирование системы и переходные процессы}

Моделирование выполнено двумя способами:

\begin{itemize}
  \item прямым программированием дискретной модели в \texttt{MATLAB} по
        уравнениям
        \[
        x_{k+1} = A_d x_k + B_d u_k,\qquad
        u_k = -K x_k + k_0 g_k,\qquad
        y_k = C_d x_k;
        \]
  \item реализацией аналогичной структуры в среде Simulink с использованием
        сумматоров, блоков усиления и задержки по одному шагу
        (\texttt{Unit Delay}), без блока \texttt{State-Space}.
\end{itemize}

На рисунке~\ref{fig:simulink_scheme} приведён вид структурной схемы
замкнутой системы в Simulink.

\begin{figure}[H]
\centering
\includegraphics[width=1\textwidth]{img/simulink_scheme.png}
\caption{Структурная схема замкнутой системы в среде Simulink}
\label{fig:simulink_scheme}
\end{figure}

Шаг моделирования для каждого варианта равен соответствующему периоду
дискретизации $T_0$, число шагов выбрано одинаковым ($N=80$), что
позволяет сравнивать длительность переходных процессов по реальному
времени $t=N T_0$.

\subsection*{Переходные процессы и управляющие воздействия}

На рисунках приведены переходные характеристики $y(k)$ и графики
управляющего воздействия $u(k)$ для обоих типов регуляторов и всех
значений периода дискретизации.

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_0_T0.1.png}
\caption{Переходная характеристика (апериодическая настройка), $T_0=0{,}1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_0_T1.0.png}
\caption{Переходная характеристика (апериодическая настройка), $T_0=1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_0_T10.0.png}
\caption{Переходная характеристика (апериодическая настройка), $T_0=10$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_30_T0.1.png}
\caption{Переходная характеристика (перерегулирование $\approx30\%$), $T_0=0{,}1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_30_T1.0.png}
\caption{Переходная характеристика (перерегулирование $\approx30\%$), $T_0=1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/resp_y_30_T10.0.png}
\caption{Переходная характеристика (перерегулирование $\approx30\%$), $T_0=10$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_0_T0.1.png}
\caption{Управляющее воздействие (апериодическая настройка), $T_0=0{,}1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_0_T1.0.png}
\caption{Управляющее воздействие (апериодическая настройка), $T_0=1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_0_T10.0.png}
\caption{Управляющее воздействие (апериодическая настройка), $T_0=10$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_30_T0.1.png}
\caption{Управляющее воздействие (перерегулирование $\approx30\%$), $T_0=0{,}1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_30_T1.0.png}
\caption{Управляющее воздействие (перерегулирование $\approx30\%$), $T_0=1$~с}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{img/ctrl_u_30_T10.0.png}
\caption{Управляющее воздействие (перерегулирование $\approx30\%$), $T_0=10$~с}
\end{figure}

% ---------- Выводы ----------
\section*{Выводы}
\addcontentsline{toc}{section}{Выводы}

В результате выполнения лабораторной работы №3:
\begin{enumerate}
  \item Для заданного непрерывного объекта управления получены дискретные
        модели при периодах дискретизации $T_0=0{,}1;\ 1;\ 10$~с. Во всех
        случаях система оказалась полностью управляемой.
  \item Синтезированы регуляторы по состоянию методом модального управления
        для двух вариантов размещения полюсов: с комплексно-сопряжённой парой
        (перерегулирование порядка $30\%$) и с вещественными полюсами
        (апериодический переходный процесс).
  \item Подбор масштабирующего коэффициента $k_0$ обеспечил нулевую
        установившуюся ошибку по ступенчатому задающему воздействию.
  \item Показано, что выбор положения полюсов позволяет управлять
        качеством переходного процесса: при комплексных полюсах достигается
        более быстрое, но колебательное регулирование, при вещественных
        полюсах переходный процесс становится апериодическим, но более
        медленным.
  \item Увеличение периода дискретизации приводит к росту коэффициентов
        обратной связи и растяжению переходных процессов во времени, что
        подтверждает необходимость разумного выбора шага дискретизации
        при проектировании дискретных систем автоматического управления.
\end{enumerate}

\appendix
\section*{Приложение А\\Листинг программы в \texttt{MATLAB}}
\addcontentsline{toc}{section}{Приложение А. Листинг программы в MATLAB}

\begin{lstlisting}[caption={Скрипт моделирования, вариант 1},label={lst:lab3}]
clc; clear; close all;

A = [ -39   -71   -23;
       23    40    13;
      -6.6 -10.7  -3.3 ];

B = [  0.228;
      -0.127;
       0.038 ];

C = [ 28000  84000  112000 ];
D = 0;

sys_c = ss(A,B,C,D);

Tvec = [0.1 1 10];

outDir = 'img';
if ~exist(outDir,'dir'), mkdir(outDir); end

fprintf('=== Дискретизация и управляемость ===\n');

for k = 1:numel(Tvec)
    T0 = Tvec(k);
    sys_d = c2d(sys_c,T0,'zoh');
    Ad = sys_d.A;
    Bd = sys_d.B;
    Cd = sys_d.C;
    Dd = sys_d.D;

    disc(k).T0 = T0;
    disc(k).Ad = Ad;
    disc(k).Bd = Bd;
    disc(k).Cd = Cd;
    disc(k).Dd = Dd;

    Wc = ctrb(Ad,Bd);
    r  = rank(Wc);

    fprintf('\nT0 = %.1f c\n',T0);
    fprintf('Ad = \n'); disp(Ad);
    fprintf('Bd = \n'); disp(Bd);
    fprintf('rank(ctrb) = %d\n',r);
end

Mp   = 0.3;
zeta = 0.36;
Ts   = 1.5;
wn   = 4/(zeta*Ts);
sigma = zeta*wn;
wd    = wn*sqrt(1 - zeta^2);

s1_30 = -sigma + 1i*wd;
s2_30 = -sigma - 1i*wd;
s3_30 = -5;

s1_0 = -2;
s2_0 = -3;
s3_0 = -4;

fprintf('\n=== Синтез регуляторов ===\n');

for k = 1:numel(Tvec)
    T0 = disc(k).T0;
    Ad = disc(k).Ad;
    Bd = disc(k).Bd;
    Cd = disc(k).Cd;

    if T0 < 2
        p30 = exp([s1_30 s2_30 s3_30]*T0);
        p0  = exp([s1_0  s2_0  s3_0 ]*T0);
    else
        p30 = [0.7+0.3j  0.7-0.3j  0.4];
        p0  = [0.7       0.5       0.3];
    end

    K30 = place(Ad,Bd,p30);
    Acl30 = Ad - Bd*K30;
    k0_30 = 1 / (Cd * ((eye(size(Ad)) - Acl30) \ Bd));

    K0 = place(Ad,Bd,p0);
    Acl0 = Ad - Bd*K0;
    k0_0 = 1 / (Cd * ((eye(size(Ad)) - Acl0) \ Bd));

    disc(k).K30   = K30;
    disc(k).K0    = K0;
    disc(k).k0_30 = k0_30;
    disc(k).k0_0  = k0_0;

    fprintf('\nT0 = %.1f c\n',T0);
    fprintf('  poles 30%% : %s\n',mat2str(p30,4));
    fprintf('  K30   = [%s]\n',num2str(K30,' %.4g'));
    fprintf('  k0_30 = %.4g\n',k0_30);
    fprintf('  poles 0%%  : %s\n',mat2str(p0,4));
    fprintf('  K0    = [%s]\n',num2str(K0,' %.4g'));
    fprintf('  k0_0  = %.4g\n',k0_0);
end

fprintf('\n=== Моделирование ===\n');

N = 80;

for k = 1:numel(Tvec)
    T0 = disc(k).T0;
    Ad = disc(k).Ad;
    Bd = disc(k).Bd;
    Cd = disc(k).Cd;

    t = (0:N-1)' * T0;
    g = ones(N,1);

    % 30% перерегулирование
    K  = disc(k).K30;
    k0 = disc(k).k0_30;

    x = zeros(3,1);
    y = zeros(N,1);
    u = zeros(N,1);

    for n = 1:N
        u(n) = -K*x + k0*g(n);
        y(n) = Cd*x;
        x    = Ad*x + Bd*u(n);
    end

    figure('Visible','off');
    plot(t,y,'LineWidth',1.5); grid on;
    xlabel('t, c'); ylabel('y(k)');
    title(sprintf('Переходная (30%%), T_0 = %.1f c',T0));
    saveas(gcf,fullfile(outDir,sprintf('resp_y_30_T%.1f.png',T0)));

    figure('Visible','off');
    plot(t,u,'LineWidth',1.5); grid on;
    xlabel('t, c'); ylabel('u(k)');
    title(sprintf('Управляющее (30%%), T_0 = %.1f c',T0));
    saveas(gcf,fullfile(outDir,sprintf('ctrl_u_30_T%.1f.png',T0)));

    % апериодическая настройка
    K  = disc(k).K0;
    k0 = disc(k).k0_0;

    x = zeros(3,1);
    y = zeros(N,1);
    u = zeros(N,1);

    for n = 1:N
        u(n) = -K*x + k0*g(n);
        y(n) = Cd*x;
        x    = Ad*x + Bd*u(n);
    end

    figure('Visible','off');
    plot(t,y,'LineWidth',1.5); grid on;
    xlabel('t, c'); ylabel('y(k)');
    title(sprintf('Переходная (апер.), T_0 = %.1f c',T0));
    saveas(gcf,fullfile(outDir,sprintf('resp_y_0_T%.1f.png',T0)));

    figure('Visible','off');
    plot(t,u,'LineWidth',1.5); grid on;
    xlabel('t, c'); ylabel('u(k)');
    title(sprintf('Управляющее (апер.), T_0 = %.1f c',T0));
    saveas(gcf,fullfile(outDir,sprintf('ctrl_u_0_T%.1f.png',T0)));
end
\end{lstlisting}

\end{document}