מודל ARMAX

פרק 1

מהו מודל ARMAX?

השם ARMAX הוא קיצור של Auto-Regressive Moving Average with eXogenous input, ובעברית: מודל אוטו־רגרסיבי עם ממוצע נע וקלט חיצוני. זהו אחד המודלים הליניאריים הנפוצים ביותר בזיהוי מערכות, והוא משלב שלושה מרכיבים: חלק אוטו־רגרסיבי שמייצג את התלות בפלטים הקודמים, חלק של קלט חיצוני (אקסוגני), וחלק של ממוצע נע (Moving Average) המתאר במדויק את דינמיקת הרעש המופעל על המערכת.

ההבדל המרכזי בין ARMAX ל־ARX הוא תוספת של פולינום שלישי, $C(q^{-1})$ , שמסנן את הרעש הלבן לפני שהוא משפיע על הפלט. בעזרתו, המודל יכול לתאר רעש "צבוע" – רעש שאינו לבן ושיש לו מבנה זמני משלו. זוהי תוספת פשוטה לכאורה, אך עם השלכות עמוקות הן מבחינה תיאורטית והן מבחינה מעשית.

הרעיון המרכזי

במודל ARX מניחים שהרעש לבן ועובר דרך אותה דינמיקה כמו הקלט ( $$1/A$$ ). במודל ARMAX מאפשרים לרעש דינמיקה משלו דרך הפולינום $$C$$ , וכך מקבלים מודל מדויק יותר ואומדנים לא מוטים – גם כאשר רעש המדידה רחוק מלהיות לבן.

מתי לבחור ARMAX?

ARMAX הוא הבחירה הטבעית כאשר התוצאות של מודל ARX אינן מספקות: כאשר השאריות עדיין מציגות מבנה דינמי ברור, כאשר אחוז ההתאמה נמוך מהרצוי, או כאשר ידוע מראש שרעש המדידה אינו לבן. במצבים כאלה, פולינום ה־ $$C$$ "סופג" את המבנה של הרעש ומאפשר ל־ $$A$$ ו־ $$B$$ לתאר את הדינמיקה האמיתית של המערכת.

פרק 2

למה צריך ARMAX? המוטיבציה הסטטיסטית

כדי להבין למה ARMAX קיים, חשוב להבין את הבעיה המרכזית שמודל ARX סובל ממנה: הטיה (Bias) בנוכחות רעש צבוע. במודל ARX מניחים שהרעש הלבן $$e(k)$$ מסונן דרך הפולינום $1/A(q^{-1})$ כדי להפוך לרעש המדידה. ההנחה הזו נוחה מאוד מבחינה אלגברית, אבל לעיתים קרובות היא רחוקה מהמציאות.

דוגמה אינטואיטיבית – מד טמפרטורה רועש

דמיינו מערכת חימום שאתם מודדים עם תרמוסטור איטי. המערכת עצמה היא מסדר ראשון ( $$A$$ מסדר 1), אבל הרעש על המדידה הוא רעש בעל פילטר משלו – הרעש הסביבתי עובר דרך החיישן והאלקטרוניקה המסננת אותו. אם תנסו לתאר את המערכת עם ARX, הפולינום $$A$$ ינסה "להתפתל" כדי להסביר גם את הדינמיקה של המערכת וגם את הדינמיקה של הרעש בו זמנית. התוצאה: $$A$$ יקבל ערכים שגויים, וגם $$B$$ יזוז בהתאם.

הפתרון של ARMAX

ARMAX מפריד באופן מפורש בין שני סוגי הדינמיקה: $$A$$ ו־ $$B$$ מתארים את הקשר הקלט־פלט של המערכת, בעוד ש־ $$C$$ מתאר את האופן שבו הרעש הלבן מסונן לרעש המדידה. כך כל פולינום ממוקד במשימתו, והאומדן הופך לעקבי ולא מוטה גם בנוכחות רעש צבוע.

ARX

מתאר את הרעש דרך $$1/A$$ . גמיש מבחינה חישובית אך מוטה כאשר הרעש אינו לבן או כאשר דינמיקת הרעש שונה מדינמיקת המערכת.

ARMAX

מתאר את הרעש דרך $$C/A$$ . גמיש יותר, מתאים למצבים שבהם הרעש מסונן דרך פילטר משלו, ומספק אומדנים לא מוטים גם בתנאים מאתגרים.

BJ (Box-Jenkins)

מפריד לחלוטין בין דינמיקת המערכת לבין דינמיקת הרעש ( $$B/F$$ ו־ $$C/D$$ ). הכי גמיש אך גם הכי קשה לאמידה.

פרק 3

המשוואה המרכזית של מודל ARMAX

עבור מערכת בדידה בזמן עם קלט יחיד ופלט יחיד (SISO), מודל ARMAX סטנדרטי נכתב בצורה:

y(k)+a_1y(k-1)+\cdots+a_{n_a}y(k-n_a) = b_1u(k-n_k)+\cdots+b_{n_b}u(k-n_k-n_b+1) +e(k)+c_1e(k-1)+\cdots+c_{n_c}e(k-n_c)

טבלת הסימונים המלאה:

סימון	משמעות
$$y(k)$$	הפלט הנמדד של המערכת בזמן המדגם $$k$$
$$u(k)$$	הקלט שהוזן למערכת בזמן המדגם $$k$$
$$e(k)$$	רעש לבן בעל תוחלת אפס וסטיית תקן $\sigma$
$$a_i$$	פרמטרים המתארים את ההשפעה של הפלטים הקודמים
$$b_i$$	פרמטרים המתארים את השפעת הקלטים הקודמים
$$c_i$$	פרמטרים של הממוצע הנע – מתארים את דינמיקת הרעש
$$n_a$$	סדר הפולינום $$A$$ – מספר ערכי הפלט הקודמים במודל
$$n_b$$	סדר הפולינום $$B$$ – מספר ערכי הקלט במודל
$$n_c$$	סדר הפולינום $$C$$ – מספר ערכי הרעש הקודמים
$$n_k$$	השהיית הקלט בדגימות (Pure Time Delay)

צורת הפולינומים עם אופרטור ההשהיה

כמו במקרה של ARX, ניתן לכתוב את המודל באמצעות אופרטור ההשהיה $q^{-1}$ . זוהי הצורה הקנונית והאלגנטית ביותר:

A(q^{-1})y(k)=B(q^{-1})u(k)+C(q^{-1})e(k)

כאשר שלושת הפולינומים הם:

A(q^{-1})=1+a_1q^{-1}+a_2q^{-2}+\cdots+a_{n_a}q^{-n_a} \] \[ B(q^{-1})=b_1q^{-n_k}+b_2q^{-(n_k+1)}+\cdots+b_{n_b}q^{-(n_k+n_b-1)} \] \[ C(q^{-1})=1+c_1q^{-1}+c_2q^{-2}+\cdots+c_{n_c}q^{-n_c}

פונקציות התמסורת השקולות

ניתן לחלץ ממשוואת ARMAX שתי פונקציות תמסורת נפרדות: אחת מהקלט לפלט, והשנייה מהרעש הלבן לפלט:

y(k)=\frac{B(q^{-1})}{A(q^{-1})}u(k)+\frac{C(q^{-1})}{A(q^{-1})}e(k)

הביטוי הראשון – $$B/A$$ – הוא פונקציית התמסורת של המערכת הדטרמיניסטית, ומתאר כיצד הקלט משפיע על הפלט. הביטוי השני – $$C/A$$ – הוא פונקציית התמסורת של הרעש, ומתאר כיצד הרעש הלבן הופך לרעש הצבוע שמופיע במדידה. ARMAX משתף את הפולינום $$A$$ בשתי הפונקציות, הנחה שהיא לרוב סבירה במערכות פיזיקליות רבות שבהן הרעש מצטבר על המצב הפנימי של המערכת לפני המדידה.

שימו לב

תנאי הכרחי ליציבות וזהות (Identifiability) של מודל ARMAX הוא ששני הפולינומים $$A$$ ו־ $$C$$ יהיו יציבים – כלומר, כל השורשים שלהם בתוך מעגל היחידה. במיוחד $$C$$ חייב להיות הפיך (Invertible), אחרת חיזוי הרעש לא יהיה יציב.

פרק 4

הבנת דינמיקת הרעש בארמקס

הלב של ARMAX הוא הפולינום $$C$$ , שמייצג את ההנחה ש"הרעש המדידה הוא לא לבן, אלא תוצאה של רעש לבן שעבר דרך פילטר ידוע". זוהי הנחה מציאותית מאוד כמעט בכל מערכת אמיתית: רעש מדידה נובע מתהליכים פיזיקליים (התרמיקה של חיישן, אלקטרוניקה של מגבר, הפרעות סביבתיות), והם כמעט תמיד מציגים מבנה זמני.

מודל אקוויואלנטי לרעש המדידה

הביטוי $\eta(k)=C(q^{-1})e(k)/A(q^{-1})$ מתאר את הרעש שמופעל על הפלט. אם נסתכל על ספקטרום ההספק שלו, נראה שאינו שטוח – יש בו תדרים חזקים יותר ותדרים חלשים יותר, בדיוק כמו ברעש טבעי. ככל ש־ $$C$$ מורכב יותר, כך הרעש המתואר יכול להיות מורכב יותר.

סדר $$n_c$$	מודל הרעש האקוויואלנטי	מתי מתאים
0	רעש לבן (מקבילים ל־ARX)	כאשר הרעש באמת לבן ולא צבוע
1	MA(1) – ממוצע נע מסדר ראשון	רעש בעל מבנה זמני קצר; שכיח מאוד
2	MA(2) – ממוצע נע מסדר שני	רעש מסונן בפילטר עם רוחב פס מוגבל
3+	MA(3+) – ממוצע נע מסדר גבוה	רעש מורכב, פילטרים מיוחדים, או מעט נדיר

הקשר בין $$C$$ לבין שאריות לבנות

המבחן הקלאסי לבדיקת איכות מודל הוא ניתוח השאריות. במודל מושלם, השאריות צריכות להיות רעש לבן. כאשר ARX מתאר היטב את הקלט־פלט אבל השאריות אינן לבנות, זה סימן ברור לכך שהפולינום $$C$$ חיוני, וצריך לעבור ל־ARMAX. במילים אחרות: ARMAX מספג את המבנה של הרעש לתוך $$C$$ , ומשאיר אחריו רק את הרעש הלבן $$e(k)$$ .

דוגמה כמותית

במערכת חימום, הפעלנו זיהוי ARX(2,2,1) וקיבלנו אחוז התאמה של 78%. ניתוח השאריות הראה אוטוקורלציה ברורה בעיכוב 1, מה שמעיד על רעש צבוע. הוספת פולינום $$C$$ מסדר 1 והמעבר ל־ARMAX(2,2,1,1) העלתה את אחוז ההתאמה ל־91%, והשאריות הפכו ללבנות. זה הסיפור הקלאסי של ARMAX בפעולה.

פרק 5

משמעות הסדרים $$n_a,n_b,n_c,n_k$$

בחירת ארבעת הסדרים של מודל ARMAX היא משימה מורכבת יותר מבחירת הסדרים ב־ARX, מכיוון שיש פרמטר נוסף – $$n_c$$ – שמשפיע על ההתנהגות הסטטיסטית של האומדן. גישה נכונה היא להתחיל ב־ARX, להבין את הדינמיקה הבסיסית, ואז להוסיף את $$C$$ בהדרגה.

$$n_a$$

קובע כמה ערכי פלט קודמים משפיעים על הפלט הנוכחי. זהה לזה שב־ARX ומשקף את סדר הדינמיקה של המערכת – מספר הקטבים בפונקציית התמסורת.

$$n_b$$

קובע כמה ערכי קלט קודמים נכנסים למודל. משקף את האופן שבו הקלט מתפרס על פני זמן ואת מספר האפסים בפונקציית התמסורת מהקלט לפלט.

$$n_c$$

חדש ב־ARMAX: קובע את סדר פולינום הממוצע הנע של הרעש. לרוב 1 או 2 מספיקים. ערכים גבוהים מ־3 הם נדירים בפועל ועלולים לגרום לבעיות זהות.

$$n_k$$

מייצג השהיה טהורה (Dead Time) בין הקלט לבין תחילת ההשפעה שלו על הפלט. זהה לזה שב־ARX וחשוב במערכות פיזיקליות עם זמן הולכה.

איך בוחרים את הסדרים בפועל?

התחל מ־ARX

אמוד מודל ARX ובדוק את השאריות. אם השאריות לבנות, אין צורך לעבור ל־ARMAX. אם יש בהן מבנה ברור – זהו סימן לכך ש־ARMAX יעזור.

שמור על $$n_a, n_b, n_k$$ מ־ARX

סדרי המערכת הדטרמיניסטית לרוב נשארים זהים. אם ARX(2,2,1) עבד טוב מבחינת ההתאמה, התחל ב־ARMAX(2,2,1,1).

הגדל את $$n_c$$ בהדרגה

התחל ב־ $$n_c=1$$ ובדוק אם השאריות הפכו ללבנות. אם לא, נסה $$n_c=2$$ . רק במקרים נדירים תזדקק ל־ $n_c \geq 3$ .

השתמש בקריטריונים סטטיסטיים

AIC ו־BIC מענישים מודלים מורכבים מדי. נסה כמה שילובים ובחר את המודל עם הקריטריון הטוב ביותר. ב־MATLAB ניתן להשתמש בפונקציה ivstruc או בסריקה ידנית.

דוגמה אינטואיטיבית

במערכת בקרת מהירות של מנוע, ידוע שהדינמיקה היא מסדר שני ( $$n_a=2$$ ). הקלט הוא מתח, ויש השהיה של דגימה אחת ( $$n_k=1$$ ). תגובת המערכת לקלט מצוינת ב־ $$n_b=2$$ . ניסוי הראה שהשאריות של ARX(2,2,1) מציגות אוטוקורלציה משמעותית בעיכוב 1, ולכן בחרנו $$n_c=1$$ . סך הכול: ARMAX(2,2,1,1).

פרק 6

אמידת פרמטרים במודל ARMAX

אמידת פרמטרים ב־ARMAX מורכבת באופן משמעותי מאמידת ARX, מכיוון שהמודל אינו ליניארי בפרמטרים. הסיבה: הרעש $$e(k)$$ אינו ידוע ישירות מהמדידות, ולכן לא ניתן לבנות את וקטור הרגרסורים בצורה סטנדרטית כפי שעושים ב־ARX.

הבעיה הליניאריות לא־ליניארית

אם ננסה לכתוב את ARMAX בצורת רגרסיה ליניארית, נקבל:

y(k)=\varphi^{T}(k)\theta+e(k)

כאשר וקטור הפרמטרים הוא:

\theta=\begin{bmatrix} -a_1 & \cdots & -a_{n_a} & b_1 & \cdots & b_{n_b} & c_1 & \cdots & c_{n_c} \end{bmatrix}^{T}

ווקטור הרגרסורים יידרש להכיל את ערכי הרעש הקודמים:

\varphi(k)=\begin{bmatrix} y(k-1) & \cdots & y(k-n_a) & u(k-n_k) & \cdots & e(k-1) & \cdots & e(k-n_c) \end{bmatrix}^{T}

הבעיה: ערכי $e(k-1), e(k-2), \ldots$ אינם ידועים, מכיוון שהרעש לא נמדד ישירות. זוהי הסיבה שלא ניתן להשתמש בריבועים פחותים פשוטים, וצריך לפנות לאלגוריתמים איטרטיביים.

שיטות אמידה עיקריות

חיזוי שגיאה (PEM)

שיטת Prediction Error Method היא הסטנדרט. המודל מחזה את הפלט, והשגיאה נמזערת באמצעות אופטימיזציה איטרטיבית – לרוב Gauss-Newton או Levenberg-Marquardt.

ריבועים פחותים מורחבים

Extended Least Squares (ELS) מעריך את $$e(k)$$ איטרטיבית משאריות הצעד הקודם, ומשתמש בערכים אלו ברגרסיה ליניארית בכל שלב. פשוט יחסית אך עלול להתבדר.

נראות מקסימלית

Maximum Likelihood הוא הגישה הסטטיסטית האידיאלית. כאשר הרעש גאוסיאני, ML שווה ערך ל־PEM ומספק תכונות סטטיסטיות אופטימליות.

אלגוריתם PEM בפירוט

PEM ממזער את הסכום הריבועי של שגיאות החיזוי. שגיאת החיזוי מוגדרת כך:

\varepsilon(k,\theta)=\frac{A(q^{-1},\theta)}{C(q^{-1},\theta)}y(k)-\frac{B(q^{-1},\theta)}{C(q^{-1},\theta)}u(k)

פונקציית המחיר היא:

V(\theta)=\frac{1}{N}\sum_{k=1}^{N}\varepsilon^{2}(k,\theta)

האומדן $\hat{\theta}$ הוא הערך שמינימלי את $V(\theta)$ . המינימיזציה מתבצעת איטרטיבית: בכל איטרציה מחשבים גרדיאנט וצעד לכיוון השיפור, עד התכנסות. בניגוד ל־ARX, אין פתרון אנליטי סגור, ולכן יש סיכוי להיתקע במינימום מקומי או להתבדר.

אתגרי האמידה האיטרטיבית

אופטימיזציה איטרטיבית רגישה לתנאים התחלתיים. ערך התחלתי גרוע עלול להוביל לאי־התכנסות, לאיטיות או למינימום מקומי. שיטה מקובלת היא אתחול חכם: להריץ קודם ARX עם אותם $$n_a, n_b, n_k$$ , ולהשתמש באומדנים האלו כנקודת התחלה ל־ARMAX, עם $$c_i$$ מאותחלים לאפסים. ב־MATLAB, פונקציית armax עושה זאת אוטומטית בעזרת InitOption.

תכונות סטטיסטיות של האומד

בהנחה שהרעש $$e(k)$$ הוא רעש לבן גאוסיאני עם תוחלת אפס, ובתנאי שהמודל בנוי נכון, האומד $\hat{\theta}$ של PEM נהנה מתכונות סטטיסטיות חזקות:

תכונה	פירוט
עקביות	כאשר $N \to \infty$ , האומד מתכנס לערך האמיתי בהסתברות 1.
התפלגות אסימפטוטית גאוסיאנית	סביב הערך האמיתי, האומד מתפלג כנורמלי, מה שמאפשר חישוב רווחי סמך.
יעילות אסימפטוטית	משיג את גבול קרמר־ראו, כלומר זוהי השגיאה התיאורטית הקטנה ביותר האפשרית.
חוסר הטיה	בניגוד ל־ARX בנוכחות רעש צבוע, ARMAX מספק אומדנים לא מוטים.

פרק 7

תהליך עבודה מומלץ לזיהוי מערכת עם ARMAX

תכנון הניסוי ואיסוף נתונים

כמו ב־ARX: קלט עשיר בתדרים (PRBS לרוב מצוין), קצב דגימה מהיר פי 5–10 מהדינמיקה המהירה, ולפחות כמה אלפי דגימות. ככל שהרעש גדול יותר וההשפעה דקה יותר, צריך יותר נתונים.

ניקוי והכנת הנתונים

הסר DC, טפל בערכים חסרים, ודא סנכרון בין קלט לפלט. חלק לאימון (66%) ובדיקה (33%). בדוק שאין רוויה (saturation) של החיישנים – זה ישבש את אמידת $$C$$ .

אמוד ARX קודם

ARX הוא תמיד נקודת המוצא. בנוסף לכך שהוא מהיר וזול חישובית, הוא מספק אומדן ראשוני של $$n_a, n_b, n_k$$ ושל הפרמטרים עצמם לאתחול ARMAX.

בדוק את שאריות ARX

הצג גרף אוטוקורלציה של השאריות. אם הן לבנות – ARX מספיק, אין צורך ב־ARMAX. אם יש מבנה ברור – המשך ל־ARMAX.

אמוד ARMAX איטרטיבית

התחל ב־ $$n_c=1$$ והגדל בהדרגה. בכל שלב, בדוק שאריות, יציבות הפולינומים והאם הקריטריון (AIC/BIC) משתפר. עצור כאשר התוספת אינה מצדיקה את הסיבוכיות.

אמת על נתונים חדשים

השוואה בין פלט המודל לבין פלט הבדיקה (cross-validation). נתח שאריות גם בסט הבדיקה. בדוק שהמודל יציב – כל הקטבים בתוך מעגל היחידה.

שימוש במודל

לאחר אימות, ניתן להשתמש במודל לחיזוי, לתכנון בקר (PID, MPC, LQG), או כבסיס לסימולציה ובחינה של מצבי הפעלה שונים.

פרק 8

אימות ובדיקת איכות המודל

אימות מודל ARMAX דומה לאימות ARX, אך עם דגש מיוחד על שני נושאים: איכות תיאור הרעש (מבחני שאריות מחמירים) ויציבות המודל (קטבים של $$A$$ וגם של $$C$$ ). שני אלה הם הצד הייחודי של ARMAX והם דורשים תשומת לב מיוחדת.

התאמת פלט (Fit %)

משווים בין הפלט האמיתי של המערכת לבין הפלט שחוזה המודל. ב־MATLAB דרך compare. ערכים מעל 85% נחשבים טובים במערכות רועשות.

ניתוח שאריות מחמיר

הבדיקה החשובה ביותר ב־ARMAX. השאריות חייבות להיות לבנות (אוטוקורלציה אפסית) ובלתי תלויות בקלט (קרוס־קורלציה אפסית). שני התנאים יחד הם תעודת איכות למודל.

יציבות הפולינומים

קטבים של $$A$$ חייבים להיות בתוך מעגל היחידה (יציבות דינמית), וכך גם השורשים של $$C$$ (הפיכות הרעש). בדיקה דרך pzmap.

קריטריונים סטטיסטיים מקובלים

בעת השוואה בין מודלים בעלי סדרים שונים, הקריטריונים הבאים מסייעים לבחור את המודל ה"נכון" ביותר תוך איזון בין דיוק לפשטות:

מדד	נוסחה	שימוש
FPE (Final Prediction Error)	$\mathrm{FPE}=V\frac{1+d/N}{1-d/N}$	השוואה בין מודלים בעלי סדרים שונים
AIC (Akaike)	$\mathrm{AIC}=\log(V)+\frac{2d}{N}$	איזון דיוק־פשטות, פנליזציה לסדרים גבוהים
BIC (Bayesian)	$\mathrm{BIC}=\log(V)+\frac{d\log(N)}{N}$	פנליזציה חזקה יותר מ־AIC על מודלים גדולים
Best Fit %	$100\left(1-\frac{\lVert y-\hat{y}\rVert}{\lVert y-\bar{y}\rVert}\right)$	מדד אינטואיטיבי להתאמת הפלט החזוי

ב־ARMAX, $$d$$ הוא $$n_a+n_b+n_c$$ – סך כל הפרמטרים הנאמדים. שימו לב שמודל ARMAX מתחרה לא רק עם מודלי ARMAX אחרים, אלא גם עם ARX. אם BIC של ARX נמוך משמעותית מזה של ARMAX – כלומר ARX פשוט יותר בלי לאבד דיוק – העדיפו את ARX.

כלל מעשי

ב־ARMAX, האימות הוא לא רק "האם המודל מתאר נכון את הקלט־פלט", אלא גם "האם המודל מתאר נכון את הרעש". מודל עם 95% התאמה אבל שאריות שאינן לבנות הוא מודל פגום. הקפידו על שני הצדדים יחד.

פרק 9

דוגמה מלאה ב־MATLAB

ב־MATLAB, ה־System Identification Toolbox מספק את הפונקציה armax שמטפלת בכל המורכבות של האמידה האיטרטיבית. הדוגמה הבאה מציגה תהליך מלא: השוואה ל־ARX, בחירת סדרים אופטימלית, אמידה, ניתוח שאריות והשוואת ביצועים.

% ===== Step 1: Prepare data =====

Ts = 0.01; % Sampling time [sec]

data = iddata(y, u, Ts); % y, u are column vectors

data.InputName = 'voltage';

data.OutputName = 'speed';

% ===== Step 2: Split into train / validation =====

N = size(data, 1);

data_train = data(1:floor(0.66*N));

data_val = data(floor(0.66*N)+1:end);

% ===== Step 3: Detrend =====

data_train = detrend(data_train);

data_val = detrend(data_val);

% ===== Step 4: Baseline ARX model =====

NN_arx = struc(1:5, 1:5, 1:3);

V_arx = arxstruc(data_train, data_val, NN_arx);

na_nb_nk = selstruc(V_arx, 'aic');

model_arx = arx(data_train, na_nb_nk);

% ===== Step 5: Check ARX residuals =====

figure; resid(data_val, model_arx);

% If residuals show structure, proceed to ARMAX

% ===== Step 6: Estimate ARMAX with grid search on nc =====

na = na_nb_nk(1); nb = na_nb_nk(2); nk = na_nb_nk(3);

best_aic = Inf; best_model = [];

for nc = 1:3

m = armax(data_train, [na nb nc nk]);

if aic(m) < best_aic

best_aic = aic(m);

best_model = m;

end

model_armax = best_model;

present(model_armax);

% ===== Step 7: Validate =====

figure; compare(data_val, model_arx, model_armax);

figure; resid(data_val, model_armax);

% ===== Step 8: Check stability =====

figure; pzmap(model_armax);

figure; bode(model_armax);

% ===== Step 9: Extract polynomials =====

[A, B, C] = polydata(model_armax);

disp('A polynomial:'); disp(A);

disp('B polynomial:'); disp(B);

disp('C polynomial:'); disp(C);

אפשרויות נוספות לאמידה מתקדמת

הפונקציה armax מקבלת ארגומנטים נוספים שמאפשרים שליטה עדינה על תהליך האמידה: אתחול ידני (InitialModel), בחירת אלגוריתם אופטימיזציה (SearchMethod – Gauss-Newton, Levenberg-Marquardt, וכו'), הגבלות על מספר איטרציות, וטולרנס התכנסות. במערכות מאתגרות במיוחד, שינוי ההגדרות הללו יכול לעשות את ההבדל בין מודל מתכנס לבין מודל שלא מתכנס.

עבור מערכות מרובות כניסות או מרובות יציאות (MIMO), השימוש דומה אך מבנה הסדרים הופך למטריצות. ב־MATLAB ניתן להשתמש ב־armax גם עבור MIMO, או לעבור למודל State Space עם n4sid לטיפול אלגנטי יותר במערכות מורכבות.

פרק 10

טעויות נפוצות וכיצד להימנע מהן

בעבודה עם ARMAX, רבות מהמלכודות של ARX חוזרות על עצמן (קלט לא מעורר, קצב דגימה גרוע, סדר גבוה מדי). אך יש גם מלכודות ייחודיות ל־ARMAX שכרוכות בטבע האיטרטיבי של האמידה ובמורכבות הנוספת של פולינום $$C$$ .

שימוש ב־ARMAX כברירת מחדל ללא צורך

ARMAX מורכב יותר מ־ARX ודורש יותר נתונים לאמידה אמינה. אם השאריות של ARX לבנות, אין צורך לעבור ל־ARMAX. אל תוסיפו סיבוכיות מיותרת.

אתחול גרוע של האופטימיזציה

האמידה האיטרטיבית רגישה לתנאים התחלתיים. אתחול אקראי או נאיבי עלול להוביל למינימום מקומי. השתמשו תמיד באמידת ARX קודם כנקודת התחלה.

$$n_c$$ גבוה מדי

ערכים של $n_c \geq 3$ נדירים בפועל ועלולים לגרום לבעיות זהות (Identifiability) – שני מודלים שונים שמסבירים את הנתונים באותה איכות. התחילו ב־ $$n_c=1$$ ובדקו אם זה מספיק.

קטבים של $$C$$ מחוץ למעגל היחידה

אם $$C$$ אינו הפיך, לא ניתן לחזות את הרעש בצורה יציבה והמודל הופך לבעייתי. ב־MATLAB, יש לוודא שהאופציה Stability מופעלת כדי לאלץ פתרון יציב.

חוסר התכנסות של האלגוריתם

אם האופטימיזציה לא מתכנסת, נסו: לשפר את האתחול, להגדיל את מספר האיטרציות המקסימלי, להחליף שיטת חיפוש (Levenberg-Marquardt לרוב יציב יותר), או לחשוב מחדש על הסדרים שבחרתם.

שכחת בדיקת קרוס־קורלציה

ב־ARMAX רבים מסתפקים בבדיקת אוטוקורלציה של השאריות, ושוכחים את הקרוס־קורלציה עם הקלט. שתי הבדיקות חיוניות: אוטוקורלציה בודקת שדינמיקת הרעש נתפסה, וקרוס־קורלציה בודקת שדינמיקת המערכת נתפסה.

הסתמכות עיוורת על תוכנה

armax ב־MATLAB מחזירה תמיד תוצאה, גם אם המודל בעייתי. בדקו תמיד: האם הקטבים יציבים? האם השאריות לבנות? האם הפרמטרים סבירים פיזיקלית? אל תקבלו את הפלט בלי לאמת.

פרק 11

ARMAX מול מודלים אחרים בזיהוי מערכות

ARMAX הוא חלק ממשפחה של מודלי קלט־פלט ליניאריים, וההבדל ביניהם נעוץ בעיקר בדרך שבה הם מתארים את הרעש. הטבלה הבאה מסכמת את ההבדלים ומסייעת לבחור את המודל המתאים למשימה.

מודל	צורה כללית	מודל הרעש	מתי להשתמש
ARX	$$Ay=Bu+e$$	פשוט, $$A$$ משותף לקלט ולרעש	מודל ראשון, רעש מתון, אמידה ליניארית
ARMAX	$$Ay=Bu+Ce$$	גמיש – פולינום $$C$$ מתאר רעש צבוע	כאשר רעש המדידה צבוע ו־ARX לא מספיק
OE (Output Error)	$y=\frac{B}{F}u+e$	רעש לבן ישירות על הפלט	זיהוי דינמיקה דטרמיניסטית טהורה
BJ (Box-Jenkins)	$y=\frac{B}{F}u+\frac{C}{D}e$	הכי גמיש – דינמיקה נפרדת לחלוטין לרעש	כשרוצים לתאר במדויק גם את הרעש
State Space	$x(k+1)=Ax(k)+Bu(k),\quad y(k)=Cx(k)+Du(k)$	תלוי בניסוח הסטוכסטי	MIMO, אינטגרציה עם בקרה מודרנית

ARMAX לעומת ARX

ההבדל העיקרי, כפי שראינו, הוא הוספת פולינום $$C$$ . ARMAX מתאים כאשר רעש המדידה אינו לבן, וכאשר ARX מציג שאריות עם מבנה. המחיר: אמידה איטרטיבית מורכבת יותר ודרישת נתונים גבוהה יותר. התועלת: אומדנים לא מוטים ותיאור ריאליסטי יותר של הדינמיקה.

ARMAX לעומת BJ

ב־BJ, פולינום $$F$$ נפרד מ־ $$D$$ , מה שמאפשר דינמיקה שונה לגמרי לרעש ולמערכת. זה גמיש יותר, אבל גם קשה יותר לאמידה ודורש יותר נתונים. ARMAX מהווה איזון טוב בין הגמישות של BJ לבין הפשטות של ARX.

כלל אצבע מעשי: התחילו ב־ARX. אם השאריות אינן לבנות אך הצורה הכללית של התגובה נראית טובה – נסו ARMAX. אם נראה שהדינמיקה של המערכת והרעש שונות מהותית – עברו ל־BJ. אם אתם עובדים עם מערכת MIMO או רוצים שילוב קל עם בקרה במרחב המצב – שקלו State Space.

פרק 12

יישומים מעשיים של מודל ARMAX

ARMAX נמצא בשימוש נרחב בתעשייה ובמחקר, במיוחד במערכות שבהן רעש המדידה משמעותי או ידוע מראש שאינו לבן. הנה כמה דוגמאות מהמציאות:

בקרת תהליכים תעשייתית

זיהוי דינמיקת כורים כימיים, מטילי מים, ומערכות חימום־קירור בנוכחות רעש סביבתי. ARMAX משמש כבסיס לבקרי MPC (Model Predictive Control).

בקרת רחפנים ורובוטיקה

במערכות עפיפה ובמערכות רובוטיות, רעש מדידה מ־IMU ומחיישנים אחרים לרוב צבוע. ARMAX מאפשר לתאר זאת ולתכנן בקרי טיסה מדויקים יותר.

אנרגיה ורשתות חכמות

חיזוי ביקושים חשמליים, תיאור דינמיקת תאים סולאריים תחת רעש מדידה, ובקרה של טורבינות רוח עם רעש סביבתי משתנה.

ביו־רפואה

מודלים פרמקוקינטיים של חולים, תגובת רמות גלוקוז לאינסולין, ודינמיקת לחץ דם – כולם מערכות עם רעש מדידה ביולוגי משמעותי שמתאים ל־ARMAX.

פיננסים וכלכלה

ARMAX הוא הרחבה של מודל ARMA המוכר מאוד מהסטטיסטיקה. בכלכלה הוא משמש לחיזוי משתנים תלויים תחת השפעות מקרו־כלכליות חיצוניות.

אקוסטיקה ועיבוד אותות

זיהוי תגובת אולמות, סינון רעש סביבתי, ובידוד אותות במערכות תקשורת. הפולינום $$C$$ מתאר באופן טבעי את הצביעה של הרעש בתעלות אקוסטיות.

דוגמה: בקרת מנוע DC ברמת רעש משמעותית

במערכת בקרת מהירות של מנוע DC, הקלט הוא מתח האספקה והפלט הוא מהירות הסיבוב הנמדדת באנקודר. אם החיישן רועש (למשל בגלל איכות אנקודר נמוכה או הפרעות אלקטרומגנטיות), ARX יתן אחוזי התאמה של 75–80%, בעוד ARMAX(2,2,1,1) יעלה את ההתאמה ל־90%+ ויספק שאריות לבנות. ההבדל הזה משמעותי כשמתכננים בקר MPC עליו.

פרק 13

שאלות נפוצות

מתי כדאי לעבור מ־ARX ל־ARMAX?

כאשר ניתוח השאריות של ARX מציג אוטוקורלציה משמעותית – כלומר השאריות אינן לבנות. זהו סימן ברור לכך שהרעש צבוע, ושמודל ARX אינו מתאר אותו כראוי. אם השאריות של ARX לבנות, אין צורך לעבור ל־ARMAX.

האם ARMAX מתאים לכל מערכת לא ליניארית?

לא. ARMAX הוא מודל ליניארי, בדיוק כמו ARX. במערכות לא ליניאריות מובהקות יש לעבור להרחבות כמו NARMAX (Nonlinear ARMAX), Hammerstein-Wiener, או רשתות נוירונים.

כמה דגימות צריך לאמידת ARMAX מוצלחת?

לרוב דרושים יותר דגימות מאשר ב־ARX, מכיוון שיש פרמטרים נוספים לאמוד והאמידה איטרטיבית. כלל אצבע: לפחות פי 20–30 ממספר הפרמטרים הכולל. עבור מודל ARMAX(2,2,1,1) זה אומר לפחות 200–300 דגימות, ועדיף הרבה יותר.

למה האלגוריתם לא מתכנס?

סיבות נפוצות: אתחול גרוע (התחילו תמיד מ־ARX), סדרים לא נכונים (במיוחד $$n_c$$ גבוה מדי), נתונים מועטים, או רעש דומיננטי מדי. נסו לשנות את שיטת החיפוש או לרדת בסדרים.

מה ההבדל בין ARMAX ל־ARMA?

ARMA מתאר סדרת נתונים יחידה ללא קלט חיצוני, ומשמש בעיקר בסטטיסטיקה ובכלכלה. ARMAX מוסיף את הקלט החיצוני $$u$$ ולכן מתאים לזיהוי מערכות דינמיות עם הפעלה ידועה.

איך בודקים שהמודל יציב?

שני תנאים: כל השורשים של $$A$$ (קטבי המערכת) חייבים להיות בתוך מעגל היחידה, וכן כל השורשים של $$C$$ . ב־MATLAB, pzmap מציג גרפית את הקטבים, וניתן לבדוק עם isstable.

האם ניתן להשתמש ב־ARMAX באופן מקוון (Online)?

קיימות גרסאות אדפטיביות של ARMAX, כמו Recursive PEM (RPEM) ו־Recursive ARMAX (RARMAX). הן יקרות יותר חישובית מ־RLS של ARX, אך מאפשרות מעקב אחרי שינויים בדינמיקה ובמבנה הרעש לאורך זמן.

מהו היחס בין ARMAX לבין מסנן Kalman?

שני המודלים מטפלים ברעש בצורה אופטימלית. ARMAX הוא מודל קלט־פלט שניתן להמיר למודל מצב (State Space), ואז ניתן לחשב את מסנן Kalman האופטימלי שלו. הקשר ההדדי הזה מאפשר לעבור בין שתי הצורות לפי צרכי היישום.

פרק 14

סיכום

מודל ARMAX מהווה צעד טבעי וחשוב מעבר למודל ARX. הוא מאפשר לתאר באופן ריאליסטי יותר את ההתנהגות של מערכות בעולם האמיתי – מערכות שבהן רעש המדידה אינו לבן, אלא מסונן דרך פילטר משלו. הוספת פולינום $$C$$ פשוטה מבחינה מתמטית, אך השלכותיה עמוקות: אומדנים לא מוטים, שאריות לבנות, ויכולת תיאור טובה יותר של דינמיקות מורכבות.

המחיר על השדרוג הוא ברור: האמידה אינה ליניארית עוד ודורשת אופטימיזציה איטרטיבית. זה אומר רגישות לתנאים התחלתיים, אפשרות להיתקע במינימום מקומי, ודרישת נתונים גבוהה יותר. אך כאשר ARX מתחיל לקרטע – השאריות אינן לבנות, ההתאמה נמוכה, או הרעש הצבוע ידוע מראש – ARMAX הוא הכלי הנכון.

שמונה כללי האצבע של ARMAX

1. תמיד התחילו מ־ARX לפני ARMAX.

2. עברו ל־ARMAX רק אם השאריות של ARX אינן לבנות.

3. השתמשו באומדן ARX כאתחול ל־ARMAX.

4. התחילו ב־

\(n_c=1\)

והגדילו רק לפי הצורך.

5. בדקו יציבות של

\(A\)

ו־

\(C\)

6. ודאו שהשאריות לבנות בסט הבדיקה, לא רק בסט האימון.

7. בדקו אוטוקורלציה וקרוס־קורלציה כאחד.

8. אם הרעש שונה לחלוטין מדינמיקת המערכת – עברו ל־BJ.

השלב הבא

לאחר שליטה ב־ARMAX, מומלץ להעמיק במודלים נוספים כמו OE (Output Error) ו־BJ (Box-Jenkins). מודלים אלה מאפשרים גמישות גדולה עוד יותר בתיאור הקשר בין דינמיקת המערכת לדינמיקת הרעש, ומתאימים למצבים מתקדמים ומאתגרים בזיהוי מערכות.

חזרה למודל ARX

סימון	משמעות
$\(y(k)\)$	הפלט הנמדד של המערכת בזמן המדגם $\(k\)$
$\(u(k)\)$	הקלט שהוזן למערכת בזמן המדגם $\(k\)$
$\(e(k)\)$	רעש לבן בעל תוחלת אפס וסטיית תקן $\sigma$
$\(a_i\)$	פרמטרים המתארים את ההשפעה של הפלטים הקודמים
$\(b_i\)$	פרמטרים המתארים את השפעת הקלטים הקודמים
$\(c_i\)$	פרמטרים של הממוצע הנע – מתארים את דינמיקת הרעש
$\(n_a\)$	סדר הפולינום $\(A\)$ – מספר ערכי הפלט הקודמים במודל
$\(n_b\)$	סדר הפולינום $\(B\)$ – מספר ערכי הקלט במודל
$\(n_c\)$	סדר הפולינום $\(C\)$ – מספר ערכי הרעש הקודמים
$\(n_k\)$	השהיית הקלט בדגימות (Pure Time Delay)