hoo2
/
SystemModling


			
							%
% !TEX TS-program = xelatex
% !TEX encoding   = UTF-8 Unicode
% !TEX spellcheck = el-GR
%
% AUTH report template for english
%
% Requires compilation with pdfLaTeX or XeLaTeX
%
% authors:
%   Χρήστος Χουτουρίδης ΑΕΜ 8997
%   cchoutou@ece.auth.gr

% Options:
%
% 1) mainlang=<language>
%    Default:  english
%    Set the default language of the document which affects hyphenations,
%    localization (section, dates, etc...)
%
%    example: \documentclass[mainlang=greek]{AUThReport}
%
% 2) <language>
%    Add hyphenation and typesetting support for other languages
%    Currently supports: english, greek, german, frenc
%
%    example: \documentclass[english, greek]{AUThReport}
%
% 3) short: Requests a shorter title for the document
%    Default: no short
%
%    example: \documentclass[short]{AUThReport}
%
\documentclass[a4paper, 11pt, mainlang=greek, english]{AUThReport/AUThReport}

\CurrentDate{\today}

% Greek report document setup suggestions
%---------------------------------
% Document configuration
\AuthorName{Χρήστος Χουτουρίδης}
\AuthorAEM{8997}
\AuthorMail{cchoutou@ece.auth.gr}

%\CoAuthorName{CoAuthor Name}
%\CoAuthorAEM{AEM}
%\CoAuthorMail{CoAuthor Mail}

% \WorkGroup{Ομάδα Χ}

\DocTitle{Εργασία 1}
\DocSubTitle{Εκτίμηση Άγνωστων Παραμέτρων - Μέθοδος Ελαχίστων Τετραγώνων}

\Department{Τμήμα ΗΜΜΥ. Τομέας Ηλεκτρονικής}
\ClassName{Προσομοίωση και Μοντελοποίηση Δυναμικών Συστημάτων}

\InstructorName{Γ. Ροβιθάκης}
\InstructorMail{rovithak@auth.gr}

\CoInstructorName{Λ. Μπίκας}
\CoInstructorMail{lnmpikas@ece.auth.gr}


% Local package requirements
%---------------------------------
%\usepackage{tabularx}
%\usepackage{array}
%\usepackage{commath}

\usepackage{amsmath, amssymb, amsfonts}
\usepackage{graphicx}
\usepackage{float}

\begin{document}

% Request a title page or header
\InsertTitle

\section{Εισαγωγή}
Η παρούσα εργασία έχει ως αντικείμενο τη μελέτη ενός απλοποιημένου δυναμικού συστήματος — συγκεκριμένα ενός γραμμικοποιημένου εκκρεμούς με ροπή εισόδου — με σκοπό την κατανόηση της συμπεριφοράς του και την εκτίμηση των φυσικών παραμέτρων του μέσω της μεθόδου των ελαχίστων τετραγώνων.
Στο πρώτο μέρος της εργασίας εξετάζεται η δυναμική του συστήματος, η οποία περιγράφεται από μία διαφορική εξίσωση δευτέρας τάξης, και προσεγγίζεται τόσο θεωρητικά όσο και αριθμητικά μέσω προσομοίωσης στο περιβάλλον MATLAB.
\par
Αρχικά, το σύστημα αναδιατυπώνεται σε μορφή κατάστασης ώστε να καταστεί κατάλληλο για ανάλυση και προσομοίωση.
Στη συνέχεια, εφαρμόζεται ημιτονική διέγερση και εξετάζεται η απόκριση του συστήματος. Ιδιαίτερο ενδιαφέρον παρουσιάζει ο τρόπος με τον οποίο το σύστημα φτάνει στη μόνιμη περιοδική του συμπεριφορά, καθώς και ο ρυθμός με τον οποίο αποκλίνουν οι ταλαντώσεις από τη μεταβατική φάση.
Η ανάλυση αυτή θέτει τη βάση για τα επόμενα μέρη της εργασίας, όπου επιχειρείται η εκτίμηση των παραμέτρων του μοντέλου με τη χρήση δεδομένων εξόδου του συστήματος.

\subsection{Παραδοτέα}
Τα παραδοτέα της εργασίας αποτελούνται από:
\begin{itemize}
    \item Την παρούσα αναφορά.
    \item Τον κατάλογο \textbf{scripts/}, που περιέχει τον κώδικα της MATLAB.
    \item Το \href{https://git.hoo2.net/hoo2/SystemModling}{σύνδεσμο} με το αποθετήριο που περιέχει όλο το project με τον κώδικα της MATLAB, της αναφοράς και τα παραδοτέα.
\end{itemize}

\section{Θέμα 1 – Μοντελοποίηση και Προσομοίωση Συστήματος Εκκρεμούς}

Το σύστημα που μελετάται περιγράφεται από τη γραμμικοποιημένη διαφορική εξίσωση:
\[
    mL^2 \ddot{q}(t) + c \dot{q}(t) + mgL q(t) = u(t)
\]
όπου $q(t)$ είναι η γωνία του εκκρεμούς, $u(t)$ η ροπή εισόδου, και $m$, $L$, $c$, $g$ φυσικές σταθερές του συστήματος.

Ορίζοντας ως διάνυσμα κατάστασης:
\[
    x(t) = \begin{bmatrix} q(t) \\ \dot{q}(t) \end{bmatrix}
\]
οι εξισώσεις κατάστασης γράφονται ως:
\[
    \dot{x}(t) = A x(t) + B u(t)
\]
όπου:
\[
    A =
    \begin{bmatrix}
        0            & 1 \\
        -\frac{g}{L} & -\frac{c}{mL^2}
    \end{bmatrix},
    \quad
    B = \begin{bmatrix}  0 \\ \frac{1}{mL^2} \end{bmatrix}
\]

Η συνάρτηση μεταφοράς του συστήματος προκύπτει από την αρχική εξίσωση ως:
\[
    G(s) = \frac{Q(s)}{U(s)} = \frac{1}{mL^2 s^2 + c s + mgL}
\]

Για την προσομοίωση χρησιμοποιήθηκε ημιτονική είσοδος $u(t) = A_0 \sin(\omega t)$, με $A_0 = 4$ και $\omega = 2$.
Οι υπόλοιπες παράμετροι ήταν $m = 0.75$, $L = 1.25$, $c = 0.15$, $g = 9.81$.
Ο αρχικός χρόνος προσομοίωσης ορίστηκε σε $20$ δευτερόλεπτα, όπως ζητείται στην εκφώνηση.

\paragraph*{Σημείωση:}
Η προσομοίωση του Θέματος 1 χρησίμευσε όχι μόνο για την παρατήρηση της δυναμικής του συστήματος, αλλά και για την παραγωγή των δεδομένων που χρησιμοποιήθηκαν στο Θέμα 2 για την εκτίμηση των παραμέτρων. Η ακρίβεια των αποτελεσμάτων του Θέματος 2 επιβεβαιώνει τη σωστή ρύθμιση και υλοποίηση της προσομοίωσης.

\subsection*{Παρατήρηση Συμπεριφοράς Συστήματος}

Βλέπουμε ότι το σύστημα παρουσιάζει περιοδική απόκριση, η οποία όμως δεν σταθεροποιείται γρήγορα.
Όπως φαίνεται στο Σχήμα~\ref{fig:20s}, το πλάτος των ταλαντώσεων συνεχίζει να μεταβάλλεται ακόμη και μετά από 20 δευτερόλεπτα, γεγονός που δείχνει ότι το σύστημα βρίσκεται ακόμη σε μεταβατική κατάσταση.

Η καθυστερημένη σύγκλιση οφείλεται κυρίως:
\begin{itemize}
    \item στην πολύ μικρή τιμή του συντελεστή απόσβεσης ($c = 0.15$), και
    \item στο γεγονός ότι η συχνότητα της εισόδου ($\omega = 2$ rad/s) είναι κοντά στη φυσική συχνότητα του εκκρεμούς ($\omega_n = \sqrt{g/L} \approx 2.8$ rad/s), προκαλώντας φαινόμενο ενίσχυσης (quasi-resonance).
\end{itemize}

Για τον λόγο αυτό, επεκτείναμε τη διάρκεια προσομοίωσης στα $90$ δευτερόλεπτα.
Όπως φαίνεται στο Σχήμα~\ref{fig:90s}, το σύστημα σταθεροποιείται τελικά σε περιοδική μορφή περίπου μετά από $50$ δευτερόλεπτα.

\InsertFigure{!ht}{1}{fig:20s}{../scripts/Prob1_responce_20s.png}{
    Απόκριση του συστήματος για $t \in [0, 20]$ sec. Η μεταβατική φάση παραμένει ενεργή.
}
\InsertFigure{!ht}{1}{fig:90s}{../scripts/Prob1_responce_90s.png}{
    Απόκριση του συστήματος για $t \in [0, 90]$ sec. Το σύστημα σταθεροποιείται σε περιοδική συμπεριφορά μετά τα $50$ sec.
}

\section{Θέμα 2 – Εκτίμηση Παραμέτρων με τη Μέθοδο Ελαχίστων Τετραγώνων}

Η εξίσωση που περιγράφει τη δυναμική του εκκρεμούς είναι της μορφής:
\[
    mL^2 \ddot{q}(t) + c \dot{q}(t) + mgL q(t) = u(t)
\]
Η παραπάνω εξίσωση είναι γραμμική ως προς τις παραμέτρους, επομένως μπορεί να αναδιατυπωθεί ως:
\[
    u(t) = \theta_1 \ddot{q}(t) + \theta_2 \dot{q}(t) + \theta_3 q(t)
    \quad \text{με} \quad
    \theta = \begin{bmatrix} mL^2 \\ c \\ mgL \end{bmatrix}
\]

Στόχος του παρόντος ερωτήματος είναι η εκτίμηση του διανύσματος $\theta$ με χρήση μετρήσεων από την έξοδο του συστήματος και της εισόδου $u(t)$, εφαρμόζοντας τη μέθοδο των ελαχίστων τετραγώνων.

\subsection{Υποερώτημα (α)}

Στο πρώτο υποερώτημα θεωρείται ότι είναι διαθέσιμες μετρήσεις όλων των μεταβλητών κατάστασης, δηλαδή $q(t)$, $\dot{q}(t)$ και $\ddot{q}(t)$, καθώς και της εισόδου $u(t)$.
Δημιουργείται λοιπόν το πρόβλημα παλινδρόμησης:
\[
    u = X \theta
    \quad \text{όπου} \quad
    X = \begin{bmatrix} \ddot{q}_1 & \dot{q}_1 & q_1 \\ \ddot{q}_2 & \dot{q}_2 & q_2 \\ \vdots & \vdots & \vdots \end{bmatrix}
\]
Η λύση του προκύπτει με τον γνωστό τύπο:
\[
    \hat{\theta} = (X^T X)^{-1} X^T u
\]

Αφού εκτιμηθούν οι παράμετροι, χρησιμοποιούνται για να υπολογιστεί η επιτάχυνση $\ddot{q}_{\text{est}}(t)$ και στη συνέχεια η ανακατασκευή της απόκρισης $q_{\text{est}}(t)$ με αριθμητική ολοκλήρωση.
Η διαδικασία υλοποιείται με ένα απλό βρόχο επανάληψης (for loop), όπου για κάθε χρονικό βήμα $k$ υπολογίζονται διαδοχικά:
\[
    \dot{q}_{\text{est}}(k) = \dot{q}_{\text{est}}(k-1) + T_s \cdot \ddot{q}_{\text{est}}(k-1)
    \]
    \[
    q_{\text{est}}(k) = q_{\text{est}}(k-1) + T_s \cdot \dot{q}_{\text{est}}(k-1)
\]

Στο Σχήμα~\ref{fig:prob2a} παρουσιάζονται η πραγματική και η εκτιμώμενη γωνία, καθώς και το σφάλμα $e_q(t)$ μεταξύ τους.

\InsertFigure{!ht}{1}{fig:prob2a}{../scripts/Prob2_20s_Ts0.1.png}{
    Αποτελέσματα εκτίμησης παραμέτρων με χρήση όλων των μεταβλητών κατάστασης.
}

\paragraph*{Συμπεράσματα:}
Η εκτίμηση παρουσιάζει υψηλή ακρίβεια, με σχετικό σφάλμα μικρότερο του 6\% για όλες τις παραμέτρους.
Ο αλγόριθμος κατάφερε να ανακατασκευάσει την απόκριση με μικρό σφάλμα, επιβεβαιώνοντας τη θεωρητική εγκυρότητα της μεθόδου.

\vspace{1em}

\subsection{Υποερώτημα (β)}

Στη δεύτερη περίπτωση θεωρούμε ότι διαθέσιμες είναι μόνο οι μετρήσεις της γωνίας $q(t)$ και της εισόδου $u(t)$.
Οι παράγωγοι $\dot{q}(t)$ και $\ddot{q}(t)$ υπολογίζονται αριθμητικά από το σήμα $q(t)$ με χρήση διαφορικών τελεστών 2ης τάξης ακρίβειας (κεντρικών διαφορών):
\[
    \dot{q}(t_k) \approx \frac{q_{k+1} - q_{k-1}}{2T_s},
    \quad
    \ddot{q}(t_k) \approx \frac{q_{k+1} - 2q_k + q_{k-1}}{T_s^2}
\]

Ακολουθείται η ίδια διαδικασία παλινδρόμησης με την περίπτωση (α), όπως και η ανακατασκευή της απόκρισης. Το Σχήμα~\ref{fig:prob2b} δείχνει τα αντίστοιχα αποτελέσματα.

\InsertFigure{!ht}{1}{fig:prob2b}{../scripts/Prob2b_20s_Ts0.1.png}{
    Αποτελέσματα εκτίμησης παραμέτρων με χρήση μόνο του $q(t)$ και του $u(t)$.
}

\paragraph*{Συμπεράσματα:}
Παρά τον περιορισμό στη διαθέσιμη πληροφορία, η εκτίμηση παρέμεινε εξαιρετικά ακριβής, με μικρές αποκλίσεις από την πραγματική τιμή. Το γεγονός αυτό δείχνει την ισχυρή ταυτοποιησιμότητα του συστήματος και τη δυνατότητα εφαρμογής της μεθόδου ακόμη και με ελλιπή δεδομένα.

\vspace{1em}

\subsection{Σύγκριση των δύο περιπτώσεων}

Ο παρακάτω πίνακας συγκρίνει τις πραγματικές και εκτιμώμενες τιμές για τα δύο υποερωτήματα:

\begin{center}
    \begin{tabular}{c|c|c|c|c|c}
        Παράμετρος & Πραγματική Τιμή & Εκτίμηση 2α & Σφάλμα 2α (\%) & Εκτίμηση 2β & Σφάλμα 2β (\%) \\
        \hline
        $mL^2$     & 1.1719          & 1.1007      & $-6.07\%$      & 1.0977      &  $-6.33\%$  \\
        $c$        & 0.1500          & 0.1690      & $+12.7\%$      & 0.1569      &  $+4.6\%$  \\
        $mgL$      & 9.1969          & 8.8157      & $-4.15\%$      & 8.8002      &  $-4.31\%$  \\
    \end{tabular}
\end{center}

Παρατηρούμε ότι και στις δύο περιπτώσεις οι εκτιμήσεις είναι ακριβείς και σταθερές, με τις αποκλίσεις να κυμαίνονται σε χαμηλά επίπεδα. Μάλιστα, το γεγονός ότι τα σφάλματα είναι αντίστοιχα ακόμα και όταν χρησιμοποιείται λιγότερη πληροφορία (περίπτωση 2β), αποτελεί ένδειξη για τη σταθερότητα και αξιοπιστία της μεθόδου. Ταυτόχρονα, αναδεικνύει την ιδιαίτερη σημασία της επιλογής σήματος διέγερσης και της ποιότητας των δεδομένων.


\end{document}