4.3 الحساب الحاسوبي (computer computation) لمعادلات فضاء الحالة (state-space) للزمن المستمر (CT) 

    إن موضوع الحساب الحاسوبي (computer computation) لمعادلات التفاضل والتكامل ذات الزمن المستمر (continuous-time differential and integral equations) واسع جدًا. وهو يتضمن قضايا عديدة: التحويل إلى متقطع (discretization)، وتعقيد الخوارزمية (complexity of an algorithm)، والكفاءة (efficiency)، والاستقرار العددي (numerical stability). ينبغي أن يعطي أفضل مخطط تحويل إلى متقطع النتيجة الأدق لحجم خطوة مُعطى. ويجب أن تكون الخوارزمية الجيدة بسيطة التطوير، وتتطلب عددًا صغيرًا من العمليات (الجمع والضرب)، وألا تكون حساسة لتغيرات المعلمات. هذه الموضوعات خارج نطاق هذا النص. نناقش في هذا القسم الفكرة الأساسية فقط واستخدام دوال MATLAB.

    اعتبر معادلة فضاء الحالة للزمن المستمر (CT state-space)

    $$ \begin{array}{l} \dot {x} (t) = \mathbf {A x} (t) + \mathbf {b u} (t) \\ y (t) = \mathbf {c x} (t) + d u (t) \\ \end{array} $$

    لنفترض أننا طُلب منا حساب المخرج $ y(t) $ المتولد بالمدخل $ u(t) $ عندما يكون $ t $ ضمن $ [0, t_f] $ والحالة الابتدائية (initial state) $ x(0) $ . الخطوة الأولى هي اختيار حجم خطوة $ T $ ثم تحويل معادلة فضاء الحالة للزمن المستمر إلى معادلة فضاء الحالة للزمن المتقطع (DT state-space) كما يلي

    $$ \begin{array}{l} y (k T) = \mathbf {c} _ {d} \mathbf {x} (k T) + d _ {d} u (k T) (4.18) \\ \mathbf {x} ((k + 1) T) = \mathbf {A} _ {d} \mathbf {x} (k T) + \mathbf {b} _ {d} u (k T) (4.19) \\ \end{array} $$

    إذا استخدمنا مخطط التحويل إلى متقطع في (4.11)، فإن $ \mathbf{A}_d = \mathbf{I} + T\mathbf{A} $ ، و $ \mathbf{b}_d = T\mathbf{b} $ ، و $ \mathbf{c}_d = \mathbf{c} $ ، و $ d_d = d $ . يمكن حساب مجموعة المعادلتين تكراريًا. نستخدم $ \mathbf{x}(0) $ و $ u(0) $ لحساب $ y(0) $ من (4.18) و $ \mathbf{x}(T) $ من (4.19). ثم نستخدم $ \mathbf{x}(T) $ و $ u(T) $ المحسوبين لحساب $ y(T) $ و $ \mathbf{x}(2T) $ . وبالاستمرار إلى الأمام، يمكننا حساب $ y(kT) $ لـ $ k = 0, 1, 2, \dots $ . بعد ذلك نُجري استيفاء (interpolation) لتوليد إشارة زمن مستمر (CT signal) $ y(t) $ من متتالية الزمن المتقطع (DT sequence) $ y(kT) $ . هذا الحساب يتضمن فقط الجمع والضرب؛ ولا يواجه مشكلة الحساسية المذكورة في القسم 4.1. لذا فإن معادلات فضاء الحالة هي الأنسب للحساب الحاسوبي.

    السؤال المتبقي هو كيفية اختيار حجم خطوة $ T $ . من الواضح أن صِغر حجم الخطوة يزيد الدقة، لكنه يتطلب أيضًا حسابًا أكثر. لذلك في التطبيق العملي لن نختار $ T $ صغيرة دون داعٍ. تنص نظرية أخذ العينات (sampling theorem) على أنه إذا كانت إشارة $ y(t) $ محدودة الحزمة (bandlimited) إلى $ \omega_{\text{max}} $ وإذا اختيرت فترة أخذ العينات $ T $ بحيث تكون أصغر من $ \pi / \omega_{\text{max}} $ ، فإن $ y(t) $ يمكن استعادتها بدقة من متتالية العينات $ y(kT) $ . لتطبيق النظرية يجب أن نعرف طيف التردد (frequency spectrum) لـ $ y(t) $ ، وهو غير متاح عادة، وأن نستخدم المستوفي المثالي (ideal interpolator)، وهو غير قابل للتحقيق فيزيائيًا. لذا ليس من السهل تطبيق نظرية أخذ العينات. أسهل طريقة لاختيار $ T $ هي بالمحاولة والخطأ. نختار اعتباطيًا $ T_1 $ و $ T_2 < T_1 $ ونجري الحساب. إذا اختلفت النتيجتان اختلافًا ملحوظًا، فإن $ T_1 $ ليست صغيرة بما يكفي؛ وهل $ T_2 $ صغيرة بما يكفي يُحدد بمقارنة نتيجتها مع نتيجة باستخدام $ T $ أصغر. نكرر الحساب حتى تصبح نتائج قيمتين متتاليتين $ T_i $ غير قابلة للتمييز. إذا كانت نتائج استخدام $ T_1 $ و $ T_2 $ الأصليتين غير قابلة للتمييز، فإن $ T_1 $ صغيرة بما يكفي. ولإيجاد أكبر $ T $ مقبولة، نزيد $ T_1 $ حتى تبدأ نتيجتها بالاختلاف عن الحساب السابق.

    نستخدم مثالًا لتوضيح الإجراء. اعتبر معادلة فضاء الحالة للزمن المستمر

    $$ \begin{array}{l} \left[ \begin{array}{l} \dot {x} _ {1} (t) \\ \dot {x} _ {2} (t) \\ \dot {x} _ {3} (t) \end{array} \right] = \left[ \begin{array}{c c c} 0 & 0 & - 0. 2 \\ 0 & - 0. 7 5 & 0. 2 5 \\ 0. 5 & - 0. 5 & 0 \end{array} \right] \left[ \begin{array}{l} x _ {1} (t) \\ x _ {2} (t) \\ x _ {3} (t) \end{array} \right] + \left[ \begin{array}{l} 0. 2 \\ 0 \\ 0 \end{array} \right] u (t) \\ y (t) = \left[ \begin{array}{l l l} 0 & 0 & 1 \end{array} \right] \left[ \begin{array}{c} x _ {1} (t) \\ x _ {2} (t) \\ x _ {3} (t) \end{array} \right] + 0 \times u (t) \\ \end{array} $$

    نستخدم دالة MATLAB ‏ $ 1 $ sim، وهي اختصار للمحاكاة الخطية (linear simulation)، لحساب خرج النظام المتولد بالمدخل $ u(t) = \sin 2t $ لـ $ t \geq 0 $ ، والحالة الابتدائية $ \mathbf{x}(0) = [0.0.2 - 0.1]^\prime $ ، حيث تشير علامة الشرطة العليا إلى النقل (transpose). نكتب في نافذة التحرير ما يلي:

    $\%$  Program 4.1 (f41.m)  
 $a = [0 0 -0.2; 0 -0.75 0.25; 0.5 -0.5 0]; b = [0.2; 0; 0];$ $c = [0 0 1]; d = 0;$ $\text{dog} = \text{ss}(a, b, c, d)$ ;  
 $t1 = 0:1:50; u1 = \sin(2^* t1); x = [0; 0.2; -0.1]$ ;  
 $y1 = lsim(\text{dog}, u1, t1, x)$ ;  
 $t2 = 0:0.1:50; u2 = \sin(2^* t2)$ ;  
 $y2 = lsim(\text{dog}, u2, t2, x)$ ;  
 $t3 = 0:0.01:50; u3 = \sin(2^* t3)$ ;  
 $y3 = lsim(\text{dog}, u3, t3, x)$ ;  
plot(t1, y1, t2, y2, :, t3, y3, -, [0 50], [0 0])  
xlabel('Time (s)'),ylabel('Amplitude')

    السطران الأولان يعبران عن $ \{\mathbf{A},\mathbf{b},\mathbf{c},d\} $ بصيغة MATLAB. السطر الثالث يعرّف النظام. نسمي النظام "dog"، وهو مُعرَّف باستخدام نموذج فضاء الحالة (state-space model)، والمشار إليه بـ ss. السطر الرابع $ t1 = 0:1:50 $ يشير إلى مجال الزمن [0,50] الذي سيتم حسابه بخطوة 1 أو، بصورة مكافئة، بحجم خطوة $ T_{1} = 1 $ والمدخل الموافق $ u1 = \sin (2^{*}t1) $ . ثم نستخدم lsim لحساب المخرج $ y1 = 1sim(dog,u1,t1,x) $ . الخط المتصل في الشكل 4.1 يُولَّد بواسطة plot(t1,y1) والذي يقوم أيضًا باستيفاء خطي (linear interpolation) بربط كل زوج من النقاط المتجاورة بخط مستقيم. ثم نكرر الحساب باختيار $ T_{2} = 0.1 $ كحجم خطوة، وتُرسم النتيجة $ y2 $ في الشكل 4.1 بخط منقّط باستخدام plot(t2,y2,':'). إنها مختلفة جدًا عن الخط المتصل. لذا فإن النتيجة باستخدام $ T_{1} = 1 $ غير مقبولة. عند هذه النقطة لا نعرف ما إذا كانت النتيجة باستخدام $ T_{2} = 0.1 $ ستكون مقبولة. نختار بعد ذلك $ T_{3} = 0.01 $ ونعيد الحساب. نحن


    الشكل 4.1 (FIGURE 4.1) المخرجات المحسوبة باستخدام $ T_{1} = 1 $ (خط متصل)، وباستخدام $ T_{2} = 0.1 $ (خط منقّط)، وباستخدام $ T_{3} = 0.01 $ (خط متقطع-منقّط).

    ثم نستخدم plot (t3, y3, '-.') لرسم النتيجة في الشكل 4.1 باستخدام الخط المتقطع-المنقّط. وهي غير قابلة للتمييز عن تلك باستخدام $ T_{2} = 0.1 $ . وبذلك نستنتج أن استجابة معادلة فضاء الحالة يمكن حسابها من معادلتها المتقطعة إذا كان حجم الخطوة 0.1 أو أصغر. لاحظ أن دوال الرسم الثلاث السابقة يمكن دمجها في واحدة كما في السطر الثاني من الأسفل في برنامج 4.1. دالة الرسم المدمجة ترسم أيضًا المحور الأفقي الذي يُولَّد بواسطة plot ([0 50], [0 0]). نحفظ برنامج 4.1 في ملف m باسم f41.m. إن كتابة f41 في نافذة الأوامر ستُظهر الشكل 4.1 في نافذة رسم.

    يحتوي MATLAB على الدالتين step(a,b,c,d) و impulse(a,b,c,d). تحسب الدالة step خرج نظام، في حالة استرخاء ابتدائي (initially relaxed)، متولدًا بمدخل خطوة (step input)، وتحسب impulse خرجًا متولدًا بمدخل نبضي (impulse input). عند استخدام كلتا الدالتين، لا حاجة لتحديد حجم الخطوة ومجال الزمن المراد حسابه كما في استخدام lsim. لاحظ أن كلا الدالتين step و impulse مبنيتان على lsim. ومع ذلك، فإنهما تتطلبان اختيارًا تكيفيًا (adaptive selection) لأحجام الخطوات وتتوقفان تلقائيًا عندما لا تكاد الاستجابات تتغير.

    لا يمكن توليد نبضة (impulse) في أي حاسوب. فكيف تُولد استجابة النبضة (impulse response) إذن؟ اعتبر

    $$ \dot {\mathbf {x}} (t) = \mathbf {A} \mathbf {x} (t) + \mathbf {b} u (t) $$

    إن تكاملها من $ t = 0 $ إلى $ t = 0+ $ ، حيث $ 0+ $ أكبر من 0 بمقدار لا نهائي الصغر، مع $ u(t) = \delta(t) $ ، يعطي

    $$ \mathbf {x} (0 +) - \mathbf {x} (0) = \mathbf {A} \int_ {0} ^ {0 +} \mathbf {x} (t) d t + \mathbf {b} \int_ {0} ^ {0 +} \delta (t) d t = \mathbf {A} \cdot \mathbf {0} + \mathbf {b} \cdot 1 = \mathbf {b} $$

    حيث إن التكامل الأول يساوي صفرًا لأن $ \mathbf{x}(t) $ لا تحتوي على نبضة، والتكامل الثاني يساوي 1. إذن فإن مدخل النبضة ينقل الحالة الابتدائية من $ \mathbf{x}(0) = \mathbf{0} $ إلى $ \mathbf{x}(0+) = \mathbf{b} $ . وبناءً عليه يمكن توليد استجابة النبضة (impulse response) (وهي استجابة الحالة الصفرية (zero-state response)) بوصفها استجابة دخل صفري (zero-input response) مُثارة بالحالة الابتدائية $ \mathbf{b} $ . بالفعل، هذه هي الطريقة التي يولد بها MATLAB استجابات النبضة. ومع ذلك، فإن impulse يولّد استجابة صحيحة فقط لمعادلة فضاء الحالة عندما $ d = 0 $ أو، بصورة أدق، يتجاهل ببساطة $ d \neq 0 $ .

    تنطبق دوال MATLAB أيضًا على دوال الانتقال (transfer functions) مثل step (n, de) و impulse (n, de)، حيث تشير n و de إلى معاملات البسط والمقام لدالة انتقال (transfer function). داخليًا، يستخدم MATLAB أولًا الدالة tf2ss لتحويل دالة انتقال إلى معادلة فضاء حالة ثم يُجري الحساب. في برنامج 4.1، استخدمنا نموذج فضاء الحالة (state-space model) لتعريف النظام كـ dog=ss(a, b, c, d). ويمكننا أيضًا تعريف النظام كـ dog=tf(n, de) باستخدام نموذج دالة الانتقال (transfer-function model). تذكر أن دوال الانتقال تصف استجابات الحالة الصفرية (zero-state responses) فقط. لذلك عندما نستخدم نموذج دالة الانتقال تُفترض الحالة الابتدائية تلقائيًا صفرًا. إذا استخدمنا نموذج دالة الانتقال في برنامج 4.1 فستكون النتيجة مختلفة لأن الحالة الابتدائية مختلفة عن الصفر. وإذا كانت الحالة الابتدائية صفرًا، فلا فرق بين استخدام نموذج ss أو نموذج tf.