Lab02: Submission version

Lab02/report/report.tex

@@ -116,13 +116,13 @@ y(t) = u(t) = \theta^{\top} u(t), \quad \text{με } \theta = [m,\, b,\, k]^{\to
 \dot{\hat{\theta}}(t) = \gamma \cdot u(t) \cdot \left(y(t) - \hat{\theta}^{\top}(t) u(t)\right)
 \]
 
-Η προσομοίωση του ``πραγματικού'' συστήματος υλοποιήθηκε μέσω της μεθόδου Runge–Kutta 4ης τάξης (RK4), ενώ οι παράγωγοι υπολογίστηκαν με διαφορικό σχήμα.
-Το script επιτρέπει επιλογή μεταξύ normalized και unnormalized εκτιμητή και διαφορετικές τιμές $\gamma$.
+Η προσομοίωση του ``πραγματικού'' συστήματος υλοποιήθηκε μέσω της μεθόδου \textbf{Runge–Kutta} 4ης τάξης (RK4), ενώ οι παράγωγοι υπολογίστηκαν με διαφορικό σχήμα.
+Το script επιτρέπει επιλογή μεταξύ \textbf{normalized} και \textbf{unnormalized} εκτιμητή και διαφορετικές τιμές $\gamma$.
 
 \subsubsection{Σύγκριση Normalized vs Unnormalized}
 
-Παρακάτω παρουσιάζονται τα αποτελέσματα για $\gamma = 0.66$ και $T=30$ sec, για τις δύο περιπτώσεις.
-Οι εκτιμήσεις είναι διαφορετικές, κυρίως στις παραμέτρους $m$ και $b$ — γεγονός που υποδηλώνει πως η χρήση normalization επηρεάζει όχι μόνο τη δυναμική, αλλά και το σημείο σύγκλισης.
+Παρακάτω παρουσιάζονται τα αποτελέσματα για \boldmath $\gamma = 0.66$ και $T=30$ \unboldmath \textbf{sec}, για τις δύο περιπτώσεις.
+Οι εκτιμήσεις είναι διαφορετικές, κυρίως στις παραμέτρους $m$ και $b$, γεγονός που υποδηλώνει πως η χρήση normalization επηρεάζει όχι μόνο τη δυναμική, αλλά και το σημείο σύγκλισης.
 
 \begin{figure}[H]
     \centering
@@ -153,7 +153,7 @@ y(t) = u(t) = \theta^{\top} u(t), \quad \text{με } \theta = [m,\, b,\, k]^{\to
 
 Από το σχήμα \ref{gradient_1a} παρατηρούμε ότι το σύστημα απαιτεί σημαντικό χρονικό διάστημα για να συγκλίνει σε σταθερές τιμές.
 Ειδικά για μικρές τιμές του $\gamma$, η σύγκλιση είναι πιο αργή αλλά ομαλότερη, η τιμή του οποίου επιλέχθηκε με δοκιμαστικές εκτελέσεις.
-Παρόλα αυτά το σύστημα παρουσίαζε μεγάλες ταλαντώσεις.
+Παρόλα αυτά το σύστημα παρουσίαζε \textbf{μεγάλες ταλαντώσεις}.
 Για το λόγο αυτό επιλέχθηκε να γίνει κανονικοποίηση στη μέθοδο κλίσης (normalized gradient update).
 Η χρήση κανονικοποίησης βελτιώνει σημαντικά τη σταθερότητα του εκτιμητή και οδηγεί σε εκτιμήσεις με μικρότερη διακύμανση, ιδιαίτερα για τις παραμέτρους $m$ και $b$.
 Ωστόσο, ακόμη και με κανονικοποίηση, το σφάλμα εκτίμησης παραμένει, γεγονός που πιθανόν οφείλεται είτε στην αριθμητική προσέγγιση των παραγώγων είτε στο ότι η είσοδος δεν διεγείρει επαρκώς όλες τις δυναμικές του συστήματος.
@@ -192,7 +192,7 @@ y(t) = u(t) = \theta^{\top} u(t), \quad \text{με } \theta = [m,\, b,\, k]^{\to
 Παρατηρούμε ότι σε αυτές τις περιπτώσεις η εκτίμηση της $b$ είναι σαφώς βελτιωμένη.
 \par
 Η διαφοροποίηση αυτή υποδεικνύει ότι η επιμονή διέγερσης (persistence of excitation) που απαιτείται για την εκτίμηση όλων των παραμέτρων δεν ικανοποιείται εξίσου από κάθε τύπο εισόδου.
-Για να επιτευχθεί συνολική παρατηρησιμότητα, ενδεχομένως να απαιτείται πιο σύνθετο ή ακόμα και στοχαστικό σήμα εισόδου (π.χ. sum-of-sines ή PRBS).
+Για να επιτευχθεί συνολική παρατηρησιμότητα, ενδεχομένως να απαιτείται πιο σύνθετο ή ακόμα και στοχαστικό σήμα εισόδου.
 
 
 \subsection{1β – Εκτίμηση παραμέτρων με μέθοδο Lyapunov}
@@ -216,15 +216,17 @@ K(\hat{\theta}) = \frac{1}{2} (y(t) - \hat{y}(t))^2 = \frac{1}{2} e^2(t)
 \dot{\hat{\theta}} = -\gamma \cdot \nabla K(\hat{\theta}) = \gamma \cdot e(t) \cdot \phi(t)
 \]
 
-Αυτή η μορφή δεν απαιτεί ξεχωριστό βοηθητικό δυναμικό μοντέλο ή υπολογισμό σφαλμάτων κατάστασης, αλλά στηρίζεται αποκλειστικά σε παρατηρήσιμα σήματα και προβλέψιμη γραμμική εξίσωση ως προς τις παραμέτρους.
+\boldmath
+Αυτή η μορφή δεν απαιτεί ξεχωριστό βοηθητικό δυναμικό μοντέλο ή υπολογισμό σφαλμάτων κατάστασης, αλλά \textbf{στηρίζεται αποκλειστικά σε παρατηρήσιμα σήματα} και προβλέψιμη γραμμική εξίσωση ως προς τις παραμέτρους.
 Η μέθοδος αυτή παρουσιάζει καλύτερη σταθερότητα σε σχέση με τη μέθοδο gradient, ειδικά όταν η είσοδος δεν πληροί πλήρως τη συνθήκη επιμονής διέγερσης.
 Επιπλέον, η εισαγωγή της παραμέτρου $\gamma$ επιτρέπει τον έλεγχο της ταχύτητας προσαρμογής, χωρίς να απαιτείται normalization.
-
+\unboldmath
 \subsubsection{Υλοποίηση Εκτιμητή}
 
+\boldmath
 Ο εκτιμητής Lyapunov υλοποιήθηκε με βάση το παραπάνω θεωρητικό μοντέλο.
 Επιπλέον, ανακατασκευάστηκε η έξοδος του συστήματος $x(t)$ με βάση τις εκτιμώμενες παραμέτρους, ώστε να συγκριθεί η πραγματική έξοδος με την εκτιμώμενη $\hat{x}(t)$ και να υπολογιστεί το σφάλμα $e_x(t) = x(t) - \hat{x}(t)$.
-
+\unboldmath
 \begin{figure}[H]
     \centering
     \begin{subfigure}[t]{0.48\textwidth}
@@ -247,7 +249,7 @@ K(\hat{\theta}) = \frac{1}{2} (y(t) - \hat{y}(t))^2 = \frac{1}{2} e^2(t)
 \]
 Παρατηρείται ότι οι εκτίμηση για το $b$ συγκλίνει σε τιμές κοντά στις πραγματικές, ενώ η εκτιμήσεις του $m$ και $k$ αποκλίνουν σημαντικά, ομοίως με τη μέθοδο gradient estimation.
 \par
-Η ανακατασκευή της εξόδου $\hat{x}(t)$ παρόλα αυτά, με βάση το εκτιμώμενο μοντέλο παρουσίασε πολύ καλή συμφωνία με την πραγματική $x(t)$, με το σφάλμα θέσης $e_x(t)$ να παραμένει κάτω από $0.6$ καθ' όλη τη διάρκεια της προσομοίωσης.
+Η ανακατασκευή της εξόδου $\hat{x}(t)$ παρόλα αυτά, με βάση το εκτιμώμενο μοντέλο \textbf{παρουσίασε πολύ καλή συμφωνία με την πραγματική} $x(t)$, με το σφάλμα θέσης $e_x(t)$ να παραμένει κάτω από $0.6$ καθ' όλη τη διάρκεια της προσομοίωσης.
 Αυτό ενισχύει την αξιοπιστία της εκτίμησης και αποδεικνύει την αποτελεσματικότητα της μεθόδου σε σταθερές συνθήκες.
 
 \subsection{1γ – Επίδραση διαταραχής στη μέτρηση x(t)}
@@ -269,7 +271,7 @@ K(\hat{\theta}) = \frac{1}{2} (y(t) - \hat{y}(t))^2 = \frac{1}{2} e^2(t)
 \end{figure}
 
 Από το σχήμα \ref{fig:disturbance} παρατηρούμε ότι η προσθήκη περιορισμένης διαταραχής στη μέτρηση της θέσης επηρεάζει ελάχιστα την εκτίμηση.
-Οι τελικές τιμές των παραμέτρων συγκλίνουν σχεδόν ταυτόσημα με την καθαρή περίπτωση:
+Οι τελικές τιμές των παραμέτρων \textbf{συγκλίνουν σχεδόν ταυτόσημα με την καθαρή περίπτωση}:
 \[
 \hat{m} = 0.8432, \quad \hat{b} = 0.2195, \quad \hat{k} = 0.2590
 \]
@@ -281,7 +283,7 @@ K(\hat{\theta}) = \frac{1}{2} (y(t) - \hat{y}(t))^2 = \frac{1}{2} e^2(t)
     \item δεν επηρεάζει τις παραγώγους που χρησιμοποιούνται στον εκτιμητή.
 \end{itemize}
 
-Συνεπώς, ο εκτιμητής Lyapunov παρουσιάζει ικανοποιητική ανοσία σε ήπια διαταραχή, γεγονός που ενισχύει τη χρησιμότητά του σε ρεαλιστικά σενάρια μέτρησης.
+Συνεπώς, ο εκτιμητής Lyapunov παρουσιάζει \textbf{ικανοποιητική ανοσία} σε ήπια διαταραχή, γεγονός που ενισχύει τη χρησιμότητά του σε ρεαλιστικά σενάρια μέτρησης.
 
 \section{Θέμα 2}
 
@@ -297,7 +299,7 @@ K(\hat{\theta}) = \frac{1}{2} (y(t) - \hat{y}(t))^2 = \frac{1}{2} e^2(t)
     \item $a_1, a_2, a_3, b$ είναι άγνωστες σταθερές παραμέτρων.
 \end{itemize}
 
-Η δομή του συστήματος επιβάλλει την κατασκευή εκτιμητή παραμέτρων σε πραγματικό χρόνο.
+Η δομή του συστήματος επιβάλλει την κατασκευή εκτιμητή παραμέτρων \textbf{σε πραγματικό χρόνο}.
 Το βασικό βήμα είναι η αναγραφή της εξίσωσης σε μορφή γραμμική ως προς τις άγνωστες παραμέτρους:
 \[
 \ddot{r}(t) = \theta^\top \phi(t)
@@ -313,18 +315,18 @@ T(z) = \ln\left(\frac{1 + z}{1 - z}\right)
 Τα σφάλματα κανονικοποιούνται με χρονικά μεταβαλλόμενα ή σταθερά φράγματα, ώστε να εξασφαλίζεται η σταθερότητα και η συμβατότητα με φυσικούς περιορισμούς του συστήματος.
 
 Ο ελεγκτής βασίζεται σε δύο επίπεδα ανάδρασης.
-Πρώτα, υπολογίζεται ένα setpoint ταχύτητας:
+Πρώτα, υπολογίζεται ένα \textbf{setpoint ταχύτητας}:
 \[
 z_1(t) = \frac{r(t) - r_d(t)}{\phi(t)}, \quad \alpha(t) = -k_1 T(z_1(t))
 \]
 όπου $r_d(t)$ είναι η επιθυμητή τροχιά και $\phi(t)$ ένα χρονικά φθίνον φράγμα σφάλματος.
 
-Στη συνέχεια, υπολογίζεται η τελική είσοδος ελέγχου:
+Στη συνέχεια, υπολογίζεται η \textbf{τελική είσοδος ελέγχου}:
 \[
 z_2(t) = \frac{\dot{r}(t) - \alpha(t)}{\rho}, \quad u(t) = -k_2 T(z_2(t))
 \]
 με σταθερό φράγμα $\rho$ και ενισχυτικό κέρδος $k_2$.
-Το σχήμα αυτό επιτρέπει υψηλή ακρίβεια παρακολούθησης (καθώς $\phi_\infty \downarrow$) και ταυτόχρονα περιορίζει τις μεταβατικές ταλαντώσεις μέσω της δομής κορεσμού.
+Το σχήμα αυτό επιτρέπει \textbf{υψηλή ακρίβεια παρακολούθησης} (καθώς $\phi_\infty \downarrow$) και ταυτόχρονα περιορίζει τις μεταβατικές ταλαντώσεις μέσω της δομής κορεσμού.
 
 Η μελέτη επεκτείνεται και σε περιπτώσεις με εξωτερικές διαταραχές $d(t)$, όπου αξιολογείται η ευρωστία του εκτιμητή και του ελεγκτή.
 
@@ -339,7 +341,6 @@ r_d(t) =
     0, & t \geq 20
 \end{cases}
 \]
-
 Για την παρακολούθηση τροχιάς χρησιμοποιείται ελεγκτής ανάδρασης που περιγράφηκε παραπάνω
 
 \begin{figure}[H]
@@ -349,14 +350,16 @@ r_d(t) =
     \label{fig:2a}
 \end{figure}
 
-Από το σχήμα \ref{fig:2a} φαίνεται ότι η σύγκλιση των παραμέτρων είναι ομαλή και σχετικά ακριβής:
+Από το σχήμα \ref{fig:2a} φαίνεται ότι η \textbf{σύγκλιση} των παραμέτρων είναι \textbf{ομαλή και σχετικά ακριβής}:
 \[
 \hat{a}_1 = 1.4970,\quad \hat{a}_2 = 1.0015,\quad \hat{a}_3 = 0.7439,\quad \hat{b} = 2.0051
 \]
+\boldmath
 με τις πραγματικές τιμές να είναι $a_1 = 2.0$, $a_2 = 1.0$, $a_3 = 0.5$, $b = 2.0$.
 Παρατηρούμε ότι οι παράμετροι $a_2$ και $b$ εκτιμώνται με υψηλή ακρίβεια, ενώ οι $a_1$ και $a_3$ εκτιμούνται με χαμηλή ακρίβεια.
 \par
 Η παρακολούθηση της επιθυμητής τροχιάς παρόλα αυτά είναι επιτυχής, με το $r(t)$ να ακολουθεί το $r_d(t)$ εντός των ορίων που ορίζει το φράγμα $\phi(t)$ καθ’ όλη τη διάρκεια της προσομοίωσης.
+\unboldmath
 
 \subsection{2β – Εκτίμηση παραμέτρων και ανακατασκευή κατάστασης με μετρήσιμα r(t), ṙ(t), u(t)}
 
@@ -385,13 +388,61 @@ r_d(t) =
     \label{fig:2b}
 \end{figure}
 
-Η εκτίμηση παραμέτρων συγκλίνει σε πολύ καλές τιμές:
+Η εκτίμηση παραμέτρων \textbf{συγκλίνει σε πολύ καλές τιμές}:
 \[
 \hat{a}_1 = 1.4970,\quad \hat{a}_2 = 1.0015,\quad \hat{a}_3 = 0.7439,\quad \hat{b} = 2.0051
 \]
 Όπως φαίνεται στο σχήμα \ref{fig:2b}, ενώ η ανακατασκευή της κατάστασης $r(t)$ είναι εύστοχη, με σφάλμα $e_r(t)$ να διατηρείται κάτω από $0.05\,\text{rad}$.
 Η ελαφριά χρονική υστέρηση και υποεκτίμηση σε ορισμένα διαστήματα οφείλονται στο ότι η ανακατασκευή γίνεται με βάση εκτιμώμενες παραμέτρους και αρχικές μηδενικές συνθήκες.
-Συνολικά, η μέθοδος Lyapunov αποδείχθηκε κατάλληλη για ταυτόχρονη εκτίμηση και παρακολούθηση σε πραγματικό χρόνο με αξιόπιστα αποτελέσματα.
+Συνολικά, η μέθοδος Lyapunov αποδείχθηκε \textbf{κατάλληλη για ταυτόχρονη εκτίμηση και παρακολούθηση σε πραγματικό χρόνο} με αξιόπιστα αποτελέσματα.
 
+\subsection{2γ – Εκτίμηση παραμέτρων παρουσία εξωτερικής διαταραχής}
+
+Σε αυτό το ερώτημα, επαναλαμβάνουμε την προηγούμενη διαδικασία, αλλά αυτή τη φορά λαμβάνεται υπόψη η παρουσία εξωτερικής διαταραχής:
+\[
+d(t) = 0.15 \sin(0.5 t)
+\]
+Η διαταραχή προστίθεται στο δυναμικό μοντέλο του συστήματος, αλλά θεωρείται άγνωστη από τον εκτιμητή.
+Έτσι, η μέθοδος Lyapunov επιχειρεί την εκτίμηση των παραμέτρων $a_1$, $a_2$, $a_3$, $b$ \textbf{χωρίς να λαμβάνει υπόψη} τη συνεισφορά του $d(t)$:
+\[
+\ddot{r}(t) = \theta^\top \phi(t) + d(t)
+\]
+
+Η ανακατασκευή της κατάστασης γίνεται όπως και στο προηγούμενο ερώτημα:
+\[
+\hat{\ddot{r}}(t) = \hat{\theta}^\top \phi(\hat{r}, \hat{\dot{r}}, u)
+\]
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=0.9\textwidth]{../scripts/output/Problem2c_estimation.png}
+    \caption{Εκτίμηση παραμέτρων και ανακατασκευή $r(t)$ με εξωτερική διαταραχή}
+    \label{fig:2c}
+\end{figure}
+
+Όπως παρατηρείται στο σχήμα \ref{fig:2c}, η εισαγωγή διαταραχής οδηγεί σε περιοδικό σφάλμα στη σύγκλιση της εκτιμώμενης κατάστασης $\hat{r}(t)$.
+Η ταλάντωση στο $e_r(t)$ συνδέεται άμεσα με τη συχνότητα και πλάτος της $d(t)$.
+
+Ωστόσο, οι τελικές εκτιμήσεις παραμέτρων \textbf{παραμένουν σχεδόν αναλλοίωτες}:
+\[
+\hat{a}_1 = 1.4996,\quad \hat{a}_2 = 0.9356,\quad \hat{a}_3 = 0.7244,\quad \hat{b} = 1.9556
+\]
+γεγονός που επιβεβαιώνει ότι η μέθοδος Lyapunov είναι σχετικά ανθεκτική σε διαταραχές χαμηλής ενέργειας ή συστηματικά περιοδικού χαρακτήρα.
+
+Το σφάλμα σε $a_2$ και $a_3$ αυξάνεται ελαφρώς σε σχέση με την περίπτωση χωρίς διαταραχή, καθώς η δομή του σήματος $d(t)$ εμφανίζεται στο σφάλμα.
+Παρ’ όλα αυτά, το συνολικό σύστημα εκτίμησης παραμένει σταθερό και προβλέψιμο.
+
+\section{Συμπεράσματα}
+
+Στην παρούσα εργασία μελετήθηκαν και συγκρίθηκαν διάφορες μέθοδοι εκτίμησης παραμέτρων δυναμικών συστημάτων.
+Στο αρχικό στάδιο της εργασίας εφαρμόστηκε η μέθοδος \textbf{gradient descent} για την εκτίμηση των παραμέτρων του γραμμικού συστήματος.
+Αν και η μέθοδος είναι απλή στον υπολογισμό και απαιτεί ελάχιστους πόρους, παρατηρήθηκε ότι η απόδοσή της εξαρτάται έντονα από τη μορφή του σήματος εισόδου και την επιλογή της σταθεράς εκμάθησης $\gamma$.
+Πιο συγκεκριμένα, \textbf{η μέθοδος απέτυχε να εκτιμήσει ταυτόχρονα και με ακρίβεια όλες τις παραμέτρους του συστήματος}.
+
+Στη συνέχεια μελετήθηκε η μέθοδος \textbf{Lyapunov} η οποία έδειξε ότι παρουσιάζει καλύτερη σταθερότητα και λιγότερη εξάρτηση από την επιλογή εισόδου και αρχικών συνθηκών.
+Στο πλαίσιο γραμμικού συστήματος, η μέθοδος απέδειξε καλή ακρίβεια και σταθερότητα, ακόμη και παρουσία ήπιου θορύβου.
+Στην περίπτωση του μη γραμμικού συστήματος roll, οι εκτιμητές πραγματικού χρόνου ενσωματώθηκαν με ελεγκτές παρακολούθησης τροχιάς και αποδείχθηκαν ικανοί να παρακολουθήσουν τις επιθυμητές κινήσεις και να προσαρμοστούν σταδιακά στις άγνωστες παραμέτρους του μοντέλου.
+Παράλληλα, αξιολογήθηκε η ευρωστία των εκτιμητών υπό την επίδραση περιοδικής εξωτερικής διαταραχής, με τα αποτελέσματα να δείχνουν μικρές αποκλίσεις στις εκτιμήσεις αλλά διατήρηση της γενικής σταθερότητας και της λειτουργικότητας της μεθόδου.
+Συνολικά, η προσέγγιση Lyapunov \textbf{αποδείχθηκε κατάλληλη για διαδικασίες ταυτοποίησης μοντέλων και ελέγχου σε πραγματικό χρόνο}, με δυνατότητες εφαρμογής σε ευρύτερα μη γραμμικά συστήματα.
 
 \end{document}

Lab02/scripts/Problem1a.m

@@ -1,4 +1,7 @@
-% Gradient estimation for mass-spring-damper system (Exercise 1a)
+%
+% Problem 1a: Gradient estimation for mass-spring-damper system (Exercise 1a)
+%
+clear
 
 % True system parameters
 m_true = 1.315;

Lab02/scripts/Problem1b.m

@@ -1,7 +1,7 @@
 %
 % Problem 1b: Lyapunov-based Parameter Estimation
 %
-
+clear
 
 % True system parameters
 m_true = 1.315;

Lab02/scripts/Problem1c.m

@@ -1,6 +1,7 @@
 %
 % Problem 1c: Effect of bounded sinusoidal disturbance on measurement x(t)
 %
+clear
 
 % True system parameters
 m_true = 1.315;

Lab02/scripts/Problem2a.m

@@ -1,6 +1,7 @@
 %
 % Problem 2a: Nonlinear system roll model parameter estimation without disturbance
 %
+clear
 
 % True system parameters
 a1 = 2.0;

Lab02/scripts/Problem2b.m

@@ -1,6 +1,7 @@
 %
 % Problem 2b: Estimation of unknown parameters using Lyapunov (r, dr, u measurable)
 %
+clear
 
 % True system parameters
 a1 = 2.0;

Lab02/scripts/Problem2c.m

@@ -0,0 +1,107 @@
+%
+% Problem 2c: Estimation under external disturbance d(t)
+%
+clear
+
+% True system parameters
+a1 = 2.0;
+a2 = 1.0;
+a3 = 0.5;
+b = 2.0;
+
+% Simulation setup
+Ts = 0.001;
+T_total = 30;
+t = 0:Ts:T_total;
+N = length(t);
+
+% Reference trajectory
+r_d = zeros(1, N);
+r_d(t >= 10 & t < 20) = pi/10;
+
+% Smooth bound phi(t)
+phi0 = 1.5;
+phi_inf = 0.05;
+lambda = 0.5;
+phi = (phi0 - phi_inf) * exp(-lambda * t) + phi_inf;
+
+% Control parameters
+k1 = 1.0;
+k2 = 1.0;
+rho = 1.0;
+
+% External disturbance
+d = 0.15 * sin(0.5 * t);
+
+% Initial conditions
+r = zeros(1, N);
+dr = zeros(1, N);
+ddr = zeros(1, N);
+
+% Estimated trajectory reconstruction
+r_hat = zeros(1, N);
+dr_hat = zeros(1, N);
+dr_hat(1) = 0; r_hat(1) = 0;
+
+% Parameter estimation setup
+theta_hat = zeros(4, N);
+theta_hat(:,1) = [1; 1; 1; 1];
+gamma = 0.66;
+
+alpha = zeros(1, N);
+u = zeros(1, N);
+
+for k = 1:N-1
+    % Control law
+    z1 = (r(k) - r_d(k)) / phi(k);
+    z1 = max(min(z1, 0.999), -0.999);
+    alpha(k) = -k1 * log((1 + z1) / (1 - z1));
+
+    z2 = (dr(k) - alpha(k)) / rho;
+    z2 = max(min(z2, 0.999), -0.999);
+    u(k) = -k2 * log((1 + z2) / (1 - z2));
+
+    % System dynamics (with disturbance)
+    phi_vec = [-dr(k); -sin(r(k)); dr(k)^2 * sin(2*r(k)); u(k)];
+    ddr(k) = a1 * phi_vec(1) + a2 * phi_vec(2) + a3 * phi_vec(3) + b * phi_vec(4) + d(k);
+    dr(k+1) = dr(k) + Ts * ddr(k);
+    r(k+1) = r(k) + Ts * dr(k);
+
+    % Parameter estimation (disturbance unmodeled)
+    y = ddr(k);
+    y_hat = theta_hat(:,k)' * phi_vec;
+    e = y - y_hat;
+    theta_hat(:,k+1) = theta_hat(:,k) + Ts * gamma * e * phi_vec;
+
+    % Reconstruct r_hat from estimated theta
+    phi_hat = [-dr_hat(k); -sin(r_hat(k)); dr_hat(k)^2 * sin(2 * r_hat(k)); u(k)];
+    dd_r_hat = theta_hat(:,k)' * phi_hat;
+    dr_hat(k+1) = dr_hat(k) + Ts * dd_r_hat;
+    r_hat(k+1) = r_hat(k) + Ts * dr_hat(k);
+end
+
+fprintf('\n2c: Final estimated parameters (with disturbance):\n');
+fprintf('a1: %.4f, a2: %.4f, a3: %.4f, b: %.4f\n', theta_hat(1,end), theta_hat(2,end), theta_hat(3,end), theta_hat(4,end));
+
+% Plot
+figure('Name', 'Problem 2c - Estimation under Disturbance', 'Position', [100, 100, 1280, 860]);
+sgtitle('Problem 2c - Estimation under External Disturbance');
+
+subplot(3,1,1);
+plot(t, theta_hat', 'LineWidth', 1.4);
+legend('a_1', 'a_2', 'a_3', 'b');
+ylabel('\theta estimates'); grid on; title('Εκτιμήσεις παραμέτρων');
+
+subplot(3,1,2);
+plot(t, r, 'b', t, r_hat, '--r', 'LineWidth', 1.4);
+legend('r(t)', 'r_{hat}(t)');
+ylabel('Roll angle [rad]'); title('Πραγματικό vs εκτιμώμενο r(t)'); grid on;
+
+subplot(3,1,3);
+plot(t, r - r_hat, 'k', 'LineWidth', 1.4);
+ylabel('e_r(t)'); xlabel('Time [s]'); title('Σφάλμα θέσης: r(t) - r̂(t)'); grid on;
+
+if ~exist('output', 'dir')
+    mkdir('output');
+end
+saveas(gcf, 'output/Problem2c_estimation.png');