مدلسازی سیستم — به زبان ساده

مدلسازی سیستم‌های دینامیکی، مکانیکی و الکتریکی اولین گام در فرایند طراحی کنترل‌کننده محسوب می‌شود. در مدلسازی سیستم، تلاش می‌شود تا یک مدل ریاضی مناسب برای توصیف عملکرد یک سیستم به دست آید و سپس با توجه به آن کنترل‌کننده را طراحی کرد. در مدلسازی سیستم، مدل ریاضی مربوط به آن را می‌توان یا با استفاده از قوانینی فیزیکی (Physical Laws) و یا با استفاده از داده‌های تجربی (Experimental Data) به دست آورد. در این مطلب قصد داریم به بیان نمایش تابع انتقال سیستم و نیز نمایش فضای حالت یک سیستم بپردازیم. در هر بخش چند مثال بسیار ساده برای مدلسازی سیستم مکانیکی و الکتریکی را حل می‌کنیم و در نهایت نحوه ایجاد این مدل‌ها را در نرم‌افزار متلب نشان می‌دهیم. مهم‌ترین دستورات متلب که در این مطلب مورد استفاده قرار می‌گیرند، ss و tf هستند.

مدلسازی سیستم دینامیکی

سیستم‌های دینامیکی سیستم‌هایی هستند که با گذر زمان و بر اساس یک قانون ثابت، دچار تغییر و یا تکامل می‌شوند. برای بسیاری از سیستم‌های فیزیکی، این قانون می‌تواند بر اساس مجموعه‌ای از معادلات دیفرانسیل مرتبه اول به صورت زیر بیان شود:

$$\dot{\text x} = \frac{\text{d}x}{\text{d}t} = f(\text x(t), \text u(t), t)$$

در معادله بالا، $$x(t)$$ بردار حالت نامیده می‌شود که مجموعه‌ای از متغیرهای تشکیل دهنده دینامیک سیستم در زمان $$t$$ را در خود دارد.

فیلم آموزش سیستم های کنترل خطی در فرادرس

کلیک کنید

$$u(t)$$ بردار ورودی‌های خارجی یا بیرونی سیستم در زمان $$t$$ است و $$f$$ تابعی (شاید غیرخطی) است که مشتق زمانی (نرخ تغییرات) بردار حالت را مشخص می‌کند.

بردار حالت $$x(t_1)$$ را در هر زمان آینده می‌توان به‌ صورت دقیق تعیین کرد. برای این کار به داشتن مقدار اولیه $$x(t_0)$$ و ورودی‌های گذشته نیاز داریم و با اعمال انتگرال از $$t_0$$ تا $$t_1$$ روی معادله بالا مقدار خروجی در هر لحظه به صورت دقیق به دست می‌آید.

اگرچه متغیرهای حالت، منحصر به‌ فرد نیستند، اما باید تعداد حداقل $$n$$ متغیر وجود داشته باشد که بتوان با آن‌ها حالت سیستم را بررسی و رفتار آینده آن را پیش‌بینی کرد. $$n$$، مرتبه سیستم (System Order) نامیده می‌شود و بُعد فضای حالت را مشخص می‌کند. مرتبه سیستم، معمولاً با تعداد عناصر ذخیره‌کننده انرژی در سیستم متناظر است.

مدلسازی سیستم با نمایش فضای حالت سیستم

نمایش بلوک دیاگرامی یک سیستم با کنترل‌کننده فیدبک در شکل زیر نشان داده شده است.

در شکل بالا، مقادیر را می‌توان به صورت زیر نام‌گذاری کرد:

فیلم آموزش سیستم های کنترل خطی در فرادرس

کلیک کنید

• $$R$$: ورودی سیستم
• $$Y$$: خروجی سیستم
• $$C$$: بلوک کنترلی
• $$P$$ یا $$\tilde{P}$$: بلوک پلنت یا سیستم مورد نظر
• $$F$$: حلقه فیدبک

در این قسمت می‌خواهیم یکی از روش‌های توصیف دینامیک و مدلسازی سیستم را معرفی کنیم که به مدل فضای حالت (State Space Model) معروف است. قصد داریم با بیان چند مثال، این روش را بیان کنیم.

مثال ۱ از مدلسازی سیستم

سیستم جرم و فنر زیر را در نظر بگیرید.

معادله دیفرانسیل توصیف کننده سیستم بالا را به دست آوردید و برای مدلسازی سیستم سپس به فرم ماتریسی بازنویسی کنید.

حل:

برای این کار باید معادلات دیفرانسیل یا ODE را به مجموعه‌ای از معادلات دیفرانسیل درجه یک تبدیل کرد. توجه کنید که متغیر x در این سیستم برابر با جابه‌جایی فنر از موقعیت اولیه یا $$x=0$$ در نظر می‌شود، نه فاصله آن از دیوار؛ زیرا فنر می‌تواند هم فشرده شود و هم دچار کشش شود.

با توجه به قانون دوم نیوتون، مجموع نیرو برابر با حاصل ضرب جرم در شتاب است. بنابراین در این سیستم داریم:

$$F = m a$$

در فرمول بالا، نیروی $$F$$ برابر با مجموع نیروی فنر و اصطکاک و نیروی خارجی در نظر گرفته می‌شود. همچنین با استفاده از قانون هوک (Hooke's Law)، نیروی فنر (Spring Force) را به صورت زیر می‌توان به دست آورد:

$$spring \; force = - kx$$

با توجه به قانون استوک (Stokes' Law)، نیروی اصطکاک (Friction Force) برابر است با:

$$friction \; force = -\rho \dot { x }$$

بنابراین در نهایت معادله زیر برای این سیستم به دست می‌آید:

$$\begin{align*} { -kx - \rho \dot{x} + u = m\ddot{x} }. \end{align*}$$

می‌توانیم مقادیر $$\dot{x}$$، $$x$$ و $$\ddot{x}$$ را به یک سمت از معادله منتقل کنیم و متغیر $$u$$ در سمت دیگر از معادله باقی بماند. سپس ضرایب معادله را نرمالیزه می‌کنیم:

$$\begin{align*} \ddot{x} + \frac{\rho}{m}\dot{x} + \frac{k}{m}x = \frac{u}{m}. \end{align*}$$

فرمول بالا، یک معادله دیفرانسیل مرتبه دو را در اختیار ما قرار می‌دهد. با بازنویسی این معادله در فرم کانونی خود، می‌توان آن را به مجموعه‌ای از معادلات دیفرانسیل درجه اول تبدیل کرد. متغیر سرعت (Velocity) را برابر با $$v = \dot{x}$$ در نظر می‌گیریم و دستگاه معادلات زیر را به دست می‌آوریم:

$$\begin{align*} \left\{ \begin{array}{lrrrrr} \dot{x} = & v & \\ \dot{v} =& (- \frac{k}{m})x &+ &(- \frac{\rho}{m})v &+& (\frac{1}{m})u \end{array} \right. . \end{align*}$$

برای مشاهده واضح‌تر فرم ماتریسی معادله ۱، می‌توان آن را به صورت زیر بازنویسی کرد:

$$\begin{align*} \left\{ \begin{array}{lrrrrr} \dot{x} =& (0)x &+ &(1) v &+& (0) u \\ \dot{v} =& (- \frac{k}{m})x &+ &(- \frac{\rho}{m})v &+& (\frac{1}{m})u \end{array} \right. . \end{align*}$$

در نهایت می‌توان معادله را به فرم ماتریسی مانند زیر نوشت:

$$\begin{align*} \left(\begin{matrix} \dot{x} \\ \dot{v} \end{matrix}\right) &= \left( \begin{matrix} 0 & 1 \\ -\frac{k}{m} & -\frac{\rho}{m} \end{matrix}\right) \left(\begin{matrix} x \\ v\end{matrix}\right) + \left(\begin{matrix} 0 \\ \frac{1}{m} \end{matrix}\right)u. \end{align*}$$

مثال ۲ از مدلسازی سیستم

مدار RL فرضی مانند شکل زیر را در نظر بگیرید.

معادله دیفرانسیل متناظر با این سیستم را به دست آورید و سپس با کاهش درجه معادلات دیفرانسیل به مرتبه یک، فرم ماتریسی آن را بنویسید.

حل:

در حالت کلی، نحوه انجام این کار بسیار شبیه به مثال یک است، اما در این شرایط باید از قانون ولتاژ کیرشهف (Kirchhoff’s Voltage Law) استفاده کنیم و جمع جبری افت ولتاژ در طول مدار را برابر با صفر قرار دهیم:

$$- V_S + V_R + V_L = 0$$

افت ولتاژ در طول مقاومت را می‌توان از قانون اهم (Ohm's Law) به دست آورد:

$$V_R = RI$$

همچنین با اعمال قانون فارادی (Faraday's Law) افت ولتاژ در طول سلف را می‌توان به صورت زیر محاسبه کرد:

$$V_L = L\dot{I}$$

بنابراین در حالت کلی می‌توان نوشت:

$${-V_S + RI + L\dot{I} = 0}$$

پس داریم:

$$\begin{align*} \dot{I} = -\frac{R}{L}I + \frac{1}{L}V_S. \end{align*}$$

این معادله درجه یک، توصیف کننده سیستم مثال ۲ است. بنابراین با استفاده از فضای حالت، می‌توان مدلسازی سیستم الکتریکی و مکانیکی را انجام داد.

مدلسازی سیستم با فضای حالت n بعدی

با توجه به مثال‌های ۱ و ۲، برای سیستمی با n حالت (State) و m ورودی (Input) نیز می‌توان نوشت:

$$\begin{align*} \text{state } x = \left( \begin{matrix} x_1 \\ \vdots \\ x_n \end{matrix}\right) \in R^n, \qquad \text{input } u = \left( \begin{matrix} u_1 \\ \vdots \\ u_m \end{matrix}\right) \in R^m, \end{align*}$$

نمایش فضای حالت خطی، برابر با مجموعه‌ای از معادلات دیفرانسیل در فرم ماتریسی است. پس سیستم را می‌توان به صورت زیر در فرم ماتریسی نوشت:

$$\begin{align*} \left( \begin{matrix} \dot{x}_1 \\ \vdots \\ \dot{x}_n \end{matrix}\right) = \left( \begin{matrix} {\Huge A} \\ \\ {\text{$n \times n$}\atop\text{matrix}}\end{matrix}\right) \left(\begin{matrix} x_1 \\ \vdots \\ x_n \end{matrix}\right) + \left( \begin{matrix} {\Huge B} \\ \\ {\text{$n \times m$}\atop\text{matrix}}\end{matrix}\right) \left(\begin{matrix} u_1 \\ \vdots \\ u_m \end{matrix}\right). \end{align*}$$

معادله بالا را می‌توان به صورت زیر خلاصه کرد:

$$\dot{x} = A_{n\times n}x + B_{n \times m}u.$$

در نمایش فضای حالت، فرمول بالا را معادله دینامیک (Dynamics Equation) سیستم می‌گویند. می‌توانیم به معادله دینامیک بالا، معادله خروجی را نیز اضافه کنیم. خروجی‌های $$P$$ را به صورت بردار $$y$$ نشان می‌دهیم و داریم:

$$\begin{align*} \text{output } y &= \left( \begin{matrix} y_1 \\ \vdots \\ y_p \end{matrix}\right) \in R^p, \qquad y = Cx, \end{align*}$$

در این فرمول، $$C$$ برابر با یک ماتریس $$P \times n$$ است. به عنوان مثال، اگر فقط حالت $$x_1$$ مورد نظر باشد، می‌توانیم فقط آن را اندازه بگیریم و بنویسیم:

$$\begin{align*} y = x_1 = \left( \begin{matrix} 1 & 0 & \ldots & 0 \end{matrix} \right)\left( \begin{matrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{matrix}\right). \end{align*}$$

بنابراین در حالت کلی، نمایش فضای حالت سیستم را می‌توان به فرم زیر نوشت:

$$\begin{align*}\dot{x} &= Ax + Bu,\\ y &= Cx. \end{align*}$$

مدلسازی سیستم بر اساس نمایش فضای حالت، برای انواع مختلف سیستم‌ها مفید و کاربردی است. از مثال متداول پاندول ساده برای توضیح بیشتر کار با مدل فضای حالت استفاده می‌کنیم.

مثال ۳ از مدلسازی سیستم

یک سیستم پاندول ساده را مطابق شکل زیر در نظر بگیرید.

معادلات دیفرانسیل توصیف‌کننده سیستم را بنویسید و سپس مدل فضای حالت آن را به دست آوردید.

حل

مشابه با مثال یک، در این مثال نیز از قانون دوم نیوتون برای حرکت دورانی (چرخشی) استفاده می‌کنیم. بر اساس این قانون، مجموع گشتاور برابر با اینرسی لحظه‌ای در شتاب زاویه‌ای است. بنابراین داریم:

$$T = J \alpha$$

در فرمول بالا، گشتاور کلی یا $$T$$، برابر با مجموع گشتاور پاندول و گشتاور خارجی است. گشتاور پاندول و اینرسی لحظه‌ای هر کدام به صورت زیر به دست می‌آیند:

$$\begin{align*} \text{pendulum torque} &= \underbrace{-m {\rm g} \sin \theta}_{\text{force}} \times \underbrace{\ell}_{\text{lever arm}},\\ \text{moment of inertia } J &= m\ell^2. \end{align*}$$

پس از جایگذاری داریم:

$$ \begin{align} \label{d2_eq2}<br /> \ddot{\theta} &= - \frac{{\rm g}}{\ell}\sin \theta + \frac{1}{m\ell^2}T_{\rm e}.<br /> \end{align} $$

فرمول بالا یک معادله غیرخطی است. از ریاضیات پایه به یاد داریم که برای زوایای کوچک می‌توان از تقریب $$\sin \theta \approx \theta$$ استفاده کرد. نمایی از تقریب سینوس زاویه برای زوایای کوچک در تصویر زیر نشان داده شده است.

بنابراین داریم:

$$ \begin{align} \label{d2_eq3}<br /> \ddot{\theta} &= - \frac{{\rm g}}{\ell} \theta + \frac{1}{m\ell^2}T_{\rm e}.<br /> \end{align} $$

مطابق آن چه در مثال یک انجام دادیم، در این مثال نیز با معرفی متغیرهای $$\theta_1 = \theta$$ و $$\theta_2 = \dot{\theta}$$ داریم:

$$\begin{align*} \dot{\theta}_2 &= - \frac{{\rm g}}{\ell} \theta + \frac{1}{m\ell^2} T_{\rm e} \\ &= -\frac{{\rm g}}{\ell}\theta_1 + \frac{1}{m\ell^2}T_{\rm e}. \end{align*}$$

بنابراین نمایش فضای حالت برای یک سیستم پاندول ساده برابر است با:

$$\begin{align*} \left( \begin{matrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{matrix} \right) = \left( \begin{matrix} 0 & 1 \\ -\frac{{\rm g}}{\ell} & 0 \end{matrix} \right) \left( \begin{matrix} {\theta}_1 \\ {\theta}_2 \end{matrix} \right) + \left( \begin{matrix} 0 \\ \frac{1}{m\ell^2} \end{matrix} \right)T_{\rm e}. \end{align*}$$

البته سوالی که ممکن است پیش بیاید این است که اگر زاویه پاندول، زاویه کوچکی نباشد، آن‌گاه باید به چه صورت عمل کرد. فرض کنید زاویه پاندول در حدود $$\pi$$ باشد. در این صورت می‌توان از تقریب $$\sin \theta \approx \pi - \theta$$ استفاده کرد.

خطی سازی فضای حالت

همان طور که در مثال سوم دیدیم، برای به دست آوردن مدل فضای حالت خطی، از تقریب $$\sin x \approx x$$ حول نقطه $$x$$ استفاده کردیم. در این بخش می‌خواهیم بدانیم که آیا یک استراتژی کلی برای خطی سازی مدل فضای حالت وجود دارد یا خیر.

فیلم آموزش سیستم های کنترل خطی در فرادرس

کلیک کنید

استفاده از سری تیلور (Taylor Series) یک روش مناسب برای خطی سازی مدل فضای حالت است. بسط سری تیلور برای متغیر $$x$$ حول نقطه $$x=x_0$$، به صورت زیر نوشته می‌شود:

$$\begin{align*} f(x) &= f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2}f''(x_0)(x-x_0)^2 + \ldots \\ &\approx f(x_0) + f'(x_0)(x-x_0). \hspace{2cm} \text{(linear approximation around $x=x_0$)} \end{align*}$$

برای یک سیستم کنترل غیر خطی در حالت کلی می‌توان نوشت:

$$\begin{align*} & \dot{x} = f(x,u) &\text{(nonlinear state-space model)} \\ & x = \left( \begin{matrix} x_1 \\ \vdots \\ x_n\end{matrix}\right), \quad u = \left( \begin{matrix} u_1 \\ \vdots \\ u_m \end{matrix}\right), \quad f = \left( \begin{matrix} f_1 \\ \vdots \\ f_n \end{matrix}\right), \end{align*}$$

می‌توان بسط سری تیلور را برای چند متغیر نیز اعمال کرد. فرض کنید $$x=0$$ و $$u=0$$ یک نقطه تعادل (Equilibrium Point) سیستم باشد، یعنی $$f(0,0) = 0$$. به عبارت دیگر، زمانی که سیستم در حال استراحت باشد و هیچ نیروی کنترلی به آن اعمال نشود، سیستم به حرکت در نمی‌آید.

حول نقطه $$(x,u)=(0,0)$$، برای هر $$f_i(x, u)$$ از $$f(x,u)$$ که $$i=1,\ldots,n$$ باشد، داریم:

$$\begin{align*} f_i(x,u) = \underbrace{f_i(0,0)}_{=0} &+ \frac{\partial f_i}{\partial x_1}(0,0)x_1 + \ldots + \frac{\partial f_i}{\partial x_n}(0,0) x_n \\ &+ \frac{\partial f_i}{\partial u_1}(0,0)u_1 + \ldots + \frac{\partial f_i}{\partial u_m}(0,0) u_m. \end{align*}$$

بعد از جمع کردن $$n$$ معادله، حال مدل فضای حالت خطی شده به صورت فرمول زیر به دست می‌آید.

$$\dot{x} = Ax + Bu, \text{ where $A_{ij} = \frac{\partial f_i}{\partial x_j} \Bigg|_{x=0 \atop u=0},\,\, B_{ik} = \frac{\partial f_i}{\partial u_k} \Bigg|_{x=0 \atop u=0}$}.$$

در معادله بالا، $$A$$ در واقع ژاکوبین $$f(x,u)$$ نسبت به متغیر $$x$$ است. به طریق مشابه، $$B$$ ژاکوبین $$f(x,u)$$ نسبت به متغیر $$u$$ در نقطه تعادل سیستم است. نکته مهمی که در اینجا باید به آن توجه کنیم این است که چون از عبارات درجه بالا صرف نظر کرده‌ایم، در نتیجه سیستم خطی به دست آمده تقریبی از سیستم واقعی است و فقط برای تغییرات کوچک حول نقطهٔ تعادل معتبر است.

حال بیاید از استراتژی کلی که در این قسمت به دست آوردیم، برای محاسبه مدل فضای حالت سیستم پاندول مثال ۳ استفاده کنیم.

حل مثال ۳ با روش بسط سری تیلور

می‌خواهیم مدل فضای حالت سیستم پاندول مثال ۳ را خطی سازی کنیم. مدل فضای حالت اصلی سیستم به صورت زیر نوشته می‌شود:

$$\begin{align*} \dot{\theta}_1 &= f_1(\theta_1,\theta_2,T_{\rm e}) = \theta_2, \hspace{5cm} \text{(already linear)}\\ \dot{\theta}_2 &= f_2(\theta_1,\theta_2,T_{\rm e}) = -\frac{{\rm g}}{\ell} \sin \theta_1 + \frac{1}{m\ell^2} T_{\rm e}. \end{align*}$$

معادله $$f_2$$ را حول نقطه تعادل $$(\theta_1,\theta_2,T_{\rm e}) = (0,0,0)$$ به صورت زیر خطی سازی می‌کنیم.

$$\begin{array}{rlrlrl} \frac{\partial f_2}{\partial \theta_1} &=\, - \frac{{\rm g}}{\ell} \cos \theta_1, & \frac{\partial f_2}{\partial \theta_2} &=\, 0, \qquad \frac{\partial f_2}{\partial T_{\rm e}} &\,= \frac{1}{m\ell^2}, \\ \frac{\partial f_2}{\partial \theta_1}\Bigg|_0 &=\, - \frac{{\rm g}}{\ell}, & \frac{\partial f_2}{\partial \theta_2}\Bigg|_0 &\,= 0, \qquad \frac{\partial f_2}{\partial T_{\rm e}}\Bigg|_0 &=\, \frac{1}{m\ell^2}. \end{array}$$

بنابراین مدل فضای حالت خطی شده سیستم پاندول برابر است با:

$$\begin{align*} \dot{\theta}_1 &= \theta_2,\\ \dot{\theta}_2 &= - \frac{{\rm g}}{\ell} \theta_1 + \frac{1}{m\ell^2}T_{\rm e}, \end{align*}$$

در نهایت مدلسازی سیستم پاندول به دست می‌آید. مشاهده می‌کنیم که معادله به دست آمده با استفاده از روش بسط سری تیلور، مشابه با جواب به دست آمده به روش تقریب زاویه $$\theta$$ است.

روش عمومی خطی سازی

به این نکته توجه کنید که اگر نقطه تعادل سیستم در نقطه‌ای به غیر از $$(x, u) = (0,0)$$ قرار داشته باشد، باید با استفاده از انتقال، نقطه تعادل غیر صفر را به مبدا بازگردانیم. نحوه این انتقال به صورت زیر است:

$$\begin{align*} &\underline{x} = x-x_0, \qquad \underline{u} = u - u_0, \\ &\underline{f}(\underline{x},\underline{u}) = f(\underline{x} + x_0, \underline{u} + u_0) = f(x,u). \end{align*}$$

چون $$\dot{\underline{x}} = \underline{f}(\underline{x},\underline{u})$$، سیستمی است که نقطه تعادل آن در مبدا قرار دارد، در نتیجه می‌توان از روش بالا برای خطی سازی مدل فضای حالت استفاده کرد. بعد از محاسبه ژاکوبین، می‌توان مدل خطی شده را بر حسب متغیرهای $$x$$ و $$u$$ به دست آورد و در نهایت از تغییر متغیر $$\underline{x} = x-x_0$$ و $$\underline{u} = u - u_0$$ استفاده کرد تا یک مدل خطی شده بر حسب متغیرهای $$x$$ و $$u$$ به دست آید.

توجه کنید که برای هر سیستم خطی، باید یک نقطه تعادل در $$(x,u) = (0,0)$$ وجود داشته باشد؛ زیرا اگر تابع $$f$$ در حال حاضر خطی باشد، آن‌گاه رابطه زیر برقرار است:

$$\begin{align*} f(x,u) = Ax + Bu \xrightarrow{x = 0,~u = 0} f(0,0) = A0 + B0 = 0. \end{align*}$$

نمایش فضای حالت با عنوان نمایش حوزه زمان نیز شناخته می‌شود. این نمایش می‌تواند به سادگی برای سیستم‌های چند ورودی - چند خروجی، سیستم‌های با شرایط اولیه غیر صفر و نیز سیستم‌های غیر خطی مورد استفاده قرار گیرد. نمایش مدل سیستم بر اساس فضای حالت، به صورت گسترده در تئوری کنترل مدرن مورد استفاده قرار می‌گیرد.

پاسخ ضربه

در بخش قبلی به استفاده از مدل فضای حالت برای مدلسازی سیستم پرداختیم. مدل فضای حالت یک روش موثر برای توصیف سیستم است.

مدل پاسخ ضربه

نترل در یک چارچوب متحد است. همان طور که گفتیم، در حالت کلی مدل فضای حالت یک سیستم را می‌توان به صورت زیر نوشت:

$$\begin{align*} \dot{x} &= Ax + Bu \\ y &= Cx, \end{align*}$$

که در فرمول بالا داریم:

• $$x(t) \in R^n$$: بردار حالت‌ها در زمان t
• $$u(t) \in R^m$$: بردار ورودی‌ها در زمان t
• $$y(t) \in R^p$$: بردار خروجی‌ها در زمان t
• $$A \in R^{n \times n}$$: ماتریس دینامیک‌های سیستم
• $$B \in R^{n \times m}$$: ماتریس کنترل یا ورودی
• $$C \in R^{p \times n}$$: ماتریس خروجی یا سنسور

حال سوالی که پیش می‌آید، این است که برای مدلسازی سیستم با یک ورودی معین $$u$$، چگونه مقدار خروجی $$y$$ را به دست آوریم؟

توجه کنید که ما در این مطلب فقط سیستم‌های تک ورودی - تک خروجی (Single-Input, Single-Output) یا SISO را در نظر می‌گیریم. به عبارت دیگر، برای تمام زمان‌های $$t$$ داریم:

$$u(t),y(t) \in R$$

در مدل فضای حالت بالا، اگر $$m=p=1$$ باشد، آن‌گاه می‌توان گفت که سیستم SISO است.

برای پاسخ به سوال مطرح شده، ابتدا باید پاسخ ضربه سیستم را تعریف کنیم. همان طور که در مطالب قبلی مجله فرادرس اشاره شد، ضربه واحد (Unit Impulse) یا تابع دلتای دیراک (Dirac’s δ-function)، یک تابع از زمان است که در شرایط زیر صدق کند:

• $$\delta(t) = 0$$ برای تمام زمان‌های $$t \neq 0$$ درست باشد.
• به ازای $$a > 0$$، عبارت $$\displaystyle \int^a_{-a}\delta(t) d t = 1$$ برقرار باشد.

در تصویر زیر نمایی از تابع دلتای دیراک یا ضربه واحد را مشاهده می‌کنید.

می‌توان تابع ضربه واحد را به صورت حدی از ضربه در نظر گرفت که مساحت زیر آن برابر با واحد باشد. این مفهوم در تصویر زیر به خوبی نشان داده شده است. وقتی $$\varepsilon \to 0$$، آن‌گاه ضربه باریک‌تر می‌شود ($$\frac{1}{\varepsilon} \to + \infty$$)، اما مساحت زیر منحنی، همان $$1$$ باقی می‌ماند.

$$\varepsilon$$ و ارتفاع $$\frac{1}{\varepsilon}$$" width="271" height="165">

بیان خروجی سیستم به صورت انتگرال کانولوشن

حال ورودی $$u(t) = \delta(t-\tau)$$ را در نظر بگیرید. این تابع در واقع بیانگر ضربه واحد در زمان $$t=\tau$$ است. اگر یک سیستم خطی و نامتغیر با زمان (Time-Invariant) یا LTI مانند تصویر زیر را در نظر بگیریم که شرایط اولیه این سیستم صفر باشد، آن‌گاه می‌توان گفت که خروجی سیستم برابر با $$y(t) = h(t-\tau)$$ خواهد بود. در این حالت، تابع $$h$$ را پاسخ ضربه سیستم می‌گویند.

با در اختیار داشتن پاسخ ضربه یک سیستم، می‌توان پاسخ آن سیستم به هر ورودی دلخواه دیگر را نیز به دست آورد. در حقیقت، با توجه به خاصیت انتقالی (Sifting Property) تابع ضربه، در $$t=\tau$$ برای هر تابع $$f$$ که خوش رفتار (Well-Behaved) باشد، داریم:

$$\begin{align} \int^\infty_{-\infty} f(t)\delta(t-\tau)d t = f(\tau). \end{align} \;\;\;\;\;\;\;\; (1)$$

معادله بالا بیان می‌کند که تابع ضربه واحد قادر است مقدار تابع $$f$$ را به $$t = \tau$$ منتقل کند. بنابراین هر تابع منظم منطقی (Reasonably Regular) را می‌توان به صورت انتگرال توابع ضربه‌ واحد بیان کرد. برای محاسبه پاسخ سیستم به سایر ورودی‌های تصادفی با استفاده از پاسخ ضربه $$h$$، می‌توان با استفاده از خاصیت انتقالی بالا، سیگنال ورودی را به فرم انتگرالی نوشت:

$$\begin{align*} u(t) = \int^\infty_{-\infty}u(\tau)\delta(t-\tau)d\tau. \end{align*}$$

برای تمام سیستم‌های خطی، پاسخ کلی سیستم که در اثر دو یا چند ورودی به وجود می‌آید، برابر با جمع جبری پاسخ سیستم به هر ورودی تکی است. این اصل را اصل جمع آثار می‌گویند.

فیلم آموزش کنترل خطی با رویکرد حل مساله در فرادرس

کلیک کنید

در نتیجه، با استفاده از اصل جمع آثار (Superposition)، پاسخ یک سیستم خطی به مجموع یا انتگرال ورودی‌ها برابر با جمع یا انتگرال پاسخ سیستم به تک تک ورودی‌ها است. در نتیجه می‌توان نوشت:

$$\large \begin{align*} u(t) &= \int^\infty_{-\infty} u(\tau)\delta(t-\tau) d\tau \quad \longrightarrow \\ y(t) &= \int^\infty_{-\infty}u(\tau)\underbrace{h(t-\tau)}_{\text{response to} \atop \text{$\delta(t-\tau)$}}d\tau. \end{align*}\;\;\;\;\;(2)$$

انتگرال در معادله بالا که مقدار $$y(t)$$ را تعیین می‌کند، کانولوشن $$y$$ و $$u$$ نام دارد. بنابراین در یک سیستم LTI با شرایط اولیه صفر، مقدار خروجی برابر با کانولوشن ورودی و پاسخ ضربه سیستم است. پس داریم:

$$\begin{align*} y(t) &= u(t) \star h(t) \\ &= h(t) \star u(t) \\ &\triangleq \int^\infty_{-\infty} u(\tau) h(t-\tau)d\tau. \end{align*}$$

اما متاسفانه فرمول بالا راه راحتی برای محاسبه مقدار خروجی $$y$$ با در اختیار داشتن ورودی $$u$$ نیست. راه حل دیگری که در عمل بسیار مورد استفاده قرار می‌گیرد، استفاده از تبدیل لاپلاس (Laplace Transform) است.

مدلسازی سیستم با تبدیل لاپلاس و تابع انتقال

در این قسمت می‌خواهیم برای مدلسازی سیستم از روش تابع انتقال استفاده کنیم. می‌دانیم که تبدیل لاپلاس دو طرفه تابع $$f(t)$$ را می‌توان به صورت زیر محاسبه کرد:

$$ \begin{align*} \tag{Laplace} \label{d3_lt}<br /> F(s) = \int^\infty_{-\infty} f(\tau)e^{-s\tau}d\tau,<br /> \end{align*} $$

در فرمول بالا، $$s \in C$$ و یک متغیر مختلط است. مقادیر سیستم در حوزه زمان و نیز در حوزه فرکانس به صورت زیر هستند:

کانولوشن در حوزه زمان به صورت $$y(t) = h(t) \star u(t)$$ است، در حالی که با گرفتن تبدیل لاپلاس از طرفین معادله بالا، می‌توان نوشت:

$$\begin{align*} Y(s) &\triangleq L \{y(t)\} \\ &= L \{h(t) \star u(t) \} \\ &= L \{h(t)\} L \{u(t)\} \hspace{2cm} \text{( equivalence btwn. time and frequency domain)}\\ &\triangleq H(s) U(s). \end{align*}$$

در نتیجه مقدار خروجی در حوزه فرکانس، که معادل با مقدار کانولوشن بیان شده در قسمت قبل است، به صورت زیر نوشته می‌شود:

$$Y(s) = H(s)U(s)$$

همچنین تبدیل لاپلاس پاسخ ضربه سیستم را می‌توان به صورت زیر محاسبه کرد:

$$H(s) = \int^\infty_{-\infty}h(\tau)e^{-s\tau}d\tau$$

تابع به دست آمده در معادله بالا را تابع انتقال (Transfer Function) سیستم می‌گویند. تابع انتقال سیستم نیز یک راه مفید برای مدلسازی سیستم است.

همان طور که گفتیم سیستم‌های LTI دارای یک مشخصه بسیار مهم هستند که بر اساس آن، اگر ورودی سیستم یک سیگنال سینوسی باشد، مقدار خروجی سیستم نیز یک سیگنالی سینوسی خواهد بود که فرکانس آن با فرکانس سیگنال ورودی برابر است و فقط مقدار دامنه و فاز آن‌ها با یکدیگر متفاوت خواهد بود. اختلاف دامنه و فاز در این دو سیگنال، تابعی از فرکانس است و بر اساس پاسخ فرکانسی سیستم به دست می‌آیند.

با استفاده از تبدیل لاپلاس، می‌توان نمایش حوزه زمان یا نمایش فضای حالت را به نمایش حوزه فرکانس و یا تابع تبدیل سیستم تبدیل کرد. این نمایش بر اساس ورودی‌ها و خروجی‌های سیستم است و تابع انتقال نام دارد. یکی دیگر از مزیت‌های استفاده از تابع انتقال برای نمایش مدل یک سیستم یا مدلسازس سیستم این است که معادلات دیفرانسیل توصیف کننده سیستم را به معادلات جبری تبدیل می‌کند که انجام محاسبات و نیز آنالیز با استفاده از آن‌ها بسیار ساده‌تر است.

اگر سیستم را علّی (Causal) فرض کنیم، در آن صورت خروجی‌ها در زمان $$t$$ از ورودی‌های زمان آینده یا $$t' > t$$ تاثیر نخواهند پذیرفت. در چنین سیستمی اگر $$h(t)$$ پاسخ ضربه سیستم در زمان $$t$$ به ورودی ضربه واحد در زمان ۰ باشد، آن‌گاه مقدار تابع پاسخ ضربه برای $$t < 0$$ برابر با صفر خواهد بود. ما تمام ورودی‌های دیگر را نیز مانند تابع ضربه واحد در $$t < 0$$ برابر با صفر فرض می‌کنیم. در این حالت می‌توانیم از تبدیل لاپلاس یک طرفه (One-Sided Laplace Transforms) استفاده کنیم و بنویسیم:

$$\begin{align*} y(t) &= \int^\infty_{\color{red}{0}} u(\tau)h(t-\tau) d\tau, \\ H(s) &= \int^\infty_{\color{red}{0}} h(\tau)e^{-s\tau}d\tau. \end{align*}$$

محاسبه تابع انتقال سیستم

در این قسمت می‌خواهیم به پرسش دیگری پاسخ دهیم. با داشتن $$u(t)$$، می‌توانیم تبدیل لاپلاس $$U(s)$$ آن را توسط نرم‌افزار متلب و یا با جدول تبدیل لاپلاس به دست آوریم. اما در مورد $$h(t)$$ و $$H(s)$$ چگونه باید عمل کرد؟ سیستم LTI نشان داده شده در قسمت‌های قبل را به یاد بیاورید.

فیلم آموزش کنترل خطی با رویکرد حل مساله در فرادرس

کلیک کنید

در این سیستم، می‌توانیم از فرمول زیر استفاده کنیم و مقدار تابع انتقال سیستم و تابع ضربه آن را به دست آوریم.

$$\begin{align*} H(s) &= C(Is - A)^{-1}B, \qquad \text{(matrix inversion)} \\ h(t) &= Ce^{At}B,\,\, t \ge 0^-. \qquad \text{(matrix exponential)} \end{align*}$$

حال می‌خواهیم چند مثال برای مدلسازی سیستم از طریق محاسبه تابع انتقال سیستم را بررسی کنیم.

مثال 1

دینامیک‌های سیستم زیر را در نظر بگیرید:

$$\begin{align*} \dot{y}(t) &= -ay(t) + u(t), &\qquad \text{( $y(t)=x(t)$, full measurement)} \\ u(t) &= e^{st}. &\qquad \text{( $u(t) = 0$ for $t < 0$)} \end{align*}$$

برای مدلسازی سیستم، تابع انتقال $$H(s)$$ سیستم را محاسبه کنید.

حل:

با داشتن $$u(t) = e^{st}, t \ge 0$$، ما S را به صورت یک عدد ثابت در نظر می‌گیریم. بنابراین داریم:

$$\begin{align*} y(t) &= \int^\infty_0 h(\tau)u(t-\tau)d\tau \qquad \text{(by $u \star h = h \star u$)}\\ &= \int^\infty_0 h(\tau)e^{s(t-\tau)}d\tau \\ &= e^{st} \int^\infty_0 h(\tau)e^{-s\tau}d\tau \\ &= e^{st}H(s). \end{align*}$$

حال با جایگذاری مقدار زیر:

$$\begin{align*} \dot{y}(t) &= \frac{d}{d t}\left(H(s)e^{st}\right) = sH(s)e^{st} \end{align*}$$

در معادله $$\dot{y}(t) = -ay(t) + u(t)$$ داریم:

$$\begin{align*}\require{cancel}<br /> sH(s)\cancel{e^{st}} &= -a H(s)\cancel{e^{st}} + \cancel{e^{st}}, \qquad \forall~ s, t > 0 \\<br /> sH(s) &= -a H(s) + 1 .<br /> \end{align*} \;\;\;\;\;\; (3) $$

اکنون با محاسبه کردن $$H(s)$$ در فرمول بالا داریم:

$$\begin{align*} H(s) &= \frac{1}{s+a}, \\ \implies y(t) &= \frac{e^{st}}{s+a}. \end{align*}\;\;\;\;\;\;\;\;\; (4)$$

اگر به سمت راست تابع انتقال بالا نگاه کنیم، آن‌گاه مقدار تابع ضربه سیستم نیز با استفاده از جدول و یا گرفتن تبدیل لاپلاس معکوس به صورت زیر به دست خواهد آمد:

$$\begin{align*} h(t) &= \begin{cases} e^{-at}, & t \ge 0 \\ 0, & t < 0 \end{cases} \end{align*}$$

این مثال را با استفاده از دیاگرام زیر می‌توان نشان داد.

در نتیجه می‌توان نوشت:

$$\begin{align*} u(t) = e^{st}, \,\, t \ge 0 \qquad \xrightarrow{x(0)=0; \text{ LTI system }} \qquad y(t) = e^{st}H(s) \end{align*}$$

مدلسازی سیستم با محاسبه پاسخ فرکانسی

تا این قسمت با نحوه محاسبه $$y(t)$$ برای یک ورودی $$u(t)$$ و پاسخ ضربه مشخص $$h$$ آشنا شدیم. در این قسمت می‌خواهیم با روش مدلسازی سیستم به کمک پاسخ فرکانسی آشنا شویم.

فیلم آموزش کنترل خطی با رویکرد حل مساله در فرادرس

کلیک کنید

با داشتن ورودی $$u(t) = e^{st}$$ در مثال یک و به ازای تمام مقادیر $$S$$ می‌خواهیم نسبت زیر را محاسبه کنیم:

$$H(s) = \frac{y(t)}{u(t)}?$$

برای محاسبه $$H(s) = \dfrac{y(t)}{u(t)}$$، باید مطابق نمودار تصویر زیر مراحل را به ازای هر مقدار S که ضروری است، تکرار کنیم.

$$H(s) = \dfrac{y(t)}{u(t)}$$" width="528" height="63">

مشکل اصلی این روش در آن است که تابع $$e^{st}$$ اگر $$s > 0$$ باشد، به سرعت افزایش می‌یابد و اگر $$s < 0$$ باشد، به سرعت به سمت صفر کاهش می‌یابد. در نتیجه، سیگنال‌های ورودی باید سیگنال‌هایی پایدار (Sustained) و محدود (Bounded) باشند. به همین دلیل است که در معادله تبدیل لاپلاس که در بالا معرفی کردیم، S مجاز است که یک عدد مختلط ($$s \in C$$) باشد.

می‌توان یک عدد مختلط به فرم $$s = a+jb \in C$$ را به صورت زیر رسم کرد.

در این حالت، a برابر با قسمت حقیقی و b برابر با قسمت موهومی عدد است. همچنین می‌توان از فرم قطبی اعداد مختلط نیز استفاده کرد. در این صورت داریم:

$$\begin{align*} s &= re^{j\varphi}, \\ r &= |s| = \sqrt{a^2 + b^2}, &\qquad\text{(magnitude or modulus)} \\ \varphi &= \angle s = \tan^{-1}\left(\frac{b}{a}\right).&\qquad\text{(phase or argument)} \end{align*}$$

همچنین بر اساس فرمول اویلر داریم:

$$e^{j\varphi} = \cos \varphi + j \sin \varphi$$

سیستم LTI با شرایط اولیه صفر را در نظر بگیرید. به این سیستم یک سیگنال ورودی سینوسی $$u(t) = A\cos(\omega t)$$ با دامنه A و فرکانس $$\omega$$ رادیان بر ثانیه داده می‌شود. حال می‌توانیم مقدار خروجی $$y(t)$$ را محاسبه کنیم. می‌دانیم که تابع کسینوس مختلط در $$z \in C$$ به صورت زیر تعریف می‌شود:

$$\cos (z) \triangleq \frac{e^{iz} + e^{-iz}}2.$$

این معادله بر اساس فرمول اویلر نیز درست است. بنابراین، به طریق خطی پاسخ $$y(t)$$ برابر است با:

$$y(t) = \frac{A}{2} \Big( H(j\omega)e^{j\omega t} + H(-j\omega)e^{-j\omega t}\Big)$$

که در این فرمول داریم:

$$\large \begin{align*} H(j\omega) &= \int^\infty_0 h(\tau)e^{-j\omega\tau}d\tau, \\ H(-j\omega) &= \int^\infty_0 \underbrace{h(\tau)e^{j\omega\tau}}_{\text{complex}\atop\text{conjugate}}d\tau \\ &= \overline{H(j\omega)}. \end{align*}$$

توجه کنید که $$h(\tau) \in R$$ یک تابع با مقدار حقیقی و $$H(j\omega) \in C$$ است. می‌توان نوشت:

$$\begin{align*} H(j\omega) &= M(\omega)e^{j\varphi(\omega)} \\ H(-j\omega) &= M(\omega)e^{-j\varphi(\omega)}, \end{align*}$$

چون هر دو تابع بالا در حوزه اعداد مختلط هستند، بنابراین اگر معادله $$y(t)$$ را با داشتن دامنه و فاز $$H(j\omega)$$ و $$H(-j\omega)$$ بازنویسی کنیم، آن‌گاه داریم:

$$\large \begin{align*} y(t) &= \frac{A}{2} M(\omega)\Big( e^{j(\omega t + \varphi(\omega))} + e^{-j(\omega t + \varphi(\omega))}\Big) \\ &= A M(\omega) \cos\big(\omega t + \varphi(\omega) \big) \hspace{3cm} \text{(steady state)} \\ &= A\underbrace{M(\omega)}_{\text{amplitude}\atop\text{magnification}} \cos\big(\omega t + \underbrace{\varphi(\omega)}_{\text{phase}\atop\text{shift}}\big). \end{align*}\;\;\;\;\;\;(5)$$

معادله بالا بیان می‌کند که پاسخ حالت پایدار یک سیستم به سیگنال کسینوسی با دامنه A و فرکانس $$\omega$$، باز هم یک سیگنال کسینوسی خواهد بود، با این تفاوت که دامنه آن به مقدار $$AM(\omega)$$ مقیاس می‌شود و فرکانس سیگنال خروجی ثابت باقی می‌ماند، اما فاز آن به $$\varphi(\omega)$$ شیفت می‌یابد. در حالت کلی می‌توان گفت که پاسخ سیستم به یک سیگنال (نه فقط یک سیگنال سینوسی) همیشه بر طبق $$Y(s) = H(s)U(s)$$ محاسبه می‌شود. برای یک سیستم با شرایط اولیه غیرصفر، پاسخ سیستم به صورت مجموع پاسخ حالت گذرا (Transient Response) و پاسخ حالت پایدار (Steady State Response) توصیف می‌شود.

مدل فضای حالت در نرم‌افزار متلب

سیستم جرم و فنر و دمپر زیر را در نظر بگیرید.

دیاگرام نیروها در سیستم جرم و فنر و دمپر نیز در تصویر زیر نشان داده شده است.

برای مدلسازی سیستم مدل فضای حالت این سیستم به صورت زیر نوشته خواهد شد:

$$\dot{\mathbf{x}}=\left[\begin{array}{c}{\dot{x}} \\ {\ddot{x}}\end{array}\right]=\left[\begin{array}{cc}{0} & {1} \\ {-\frac{k}{m}} & {-\frac{b}{m}}\end{array}\right]\left[\begin{array}{c}{x} \\ {\dot{x}}\end{array}\right]+\left[\begin{array}{c}{0} \\ {\frac{1}{m}}\end{array}\right] F(t)$$

$$y=\left[\begin{array}{ll}{1} & {0}\end{array}\right]\left[\begin{array}{l}{x} \\ {\dot{x}}\end{array}\right]$$

حال می‌خواهیم نشان دهیم که مدل فضای حالت بالا را چگونه در m-file نرم‌افزار متلب وارد کنیم.

فیلم آموزش کاربرد متلب در کنترل خطی در فرادرس

کلیک کنید

ابتدا مقادیر زیر را به هر یک از المان‌های سیستم اختصاص می‌دهیم.

یک m-file جدید در متلب باز می‌کنیم و دستورات زیر را در آن وارد می‌کنیم.

1m = 1;
2k = 1;
3b = 0.2;
4F = 1;
5
6A = [0 1; -k/m -b/m];
7B = [0 1/m]';
8C = [1 0];
9D = [0];
10
11sys = ss(A,B,C,D)

sys =
 
  A = 
         x1    x2
   x1     0     1
   x2    -1  -0.2
 
  B = 
       u1
   x1   0
   x2   1
 
  C = 
       x1  x2
   y1   1   0
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

در نتیجه توانستیم مدلسازی سیستم با استفاده از فضای حالت را در نرم‌افزار متلب ایجاد کنیم.

وارد کردن تابع انتقال در متلب

تبدیل لاپلاس یک سیستم با شرایط اولیه برابر است با:

$$m s^{2} X(s)+b s X(s)+k X(s)=F(s)$$

بنابراین تابع انتقال از نیروی ورودی به جابه‌جایی به صورت زیر به دست می‌آید:

$$\frac{X(s)}{F(s)}=\frac{1}{m s^{2}+b s+k}$$

حال برای وارد کردن تابع انتقال یک سیستم در m-file متلب، باید قطعه کد زیر را در فایلی که پارامترهای سیستم در آن تعریف شده‌اند، نوشت.

فیلم آموزش کاربرد متلب در کنترل خطی در فرادرس

کلیک کنید

1s = tf('s');
2sys = 1/(m*s^2+b*s+k)

sys =
 
         1
  ---------------
  s^2 + 0.2 s + 1
 
Continuous-time transfer function.

البته می‌توان برای تعریف کردن تابع انتقال و مدلسازی سیستم در متلب و یا سیمیولینک، از روش زیر نیز استفاده کرد.

1num = [1];
2den = [m b k];
3sys = tf(num,den)

sys =
 
         1
  ---------------
  s^2 + 0.2 s + 1
 
Continuous-time transfer function.

در این مطلب به بیان روش‌‌های مختلف مدلسازی سیستم پرداختیم و در نهایت دو مدل را در نرم افزار متلب وارد کردیم. اگر این مطلب برایتان مفید بوده است، آموزش‌های زیر نیز به شما پیشنهاد می‌شوند:

• مجموعه آموزش‌های مهندسی برق
• آموزش سیستم‌های کنترل خطی
• مجموعه آموزش‌های مهندسی کنترل
• آموزش تجزیه و تحلیل سیگنال ها و سیستم ها‎
• فیدبک (Feedback) در سیستم های کنترل — مفاهیم اصلی
• تبدیل z — به زبان ساده
• مکان هندسی ریشه ها (Root Locus) در مهندسی کنترل — به زبان ساده

^^