سیستم غیر خطی — راهنمای جامع

۳۹۴۲ بازدید
آخرین به‌روزرسانی: ۲۳ اردیبهشت ۱۴۰۲
زمان مطالعه: ۸ دقیقه
سیستم غیر خطی — راهنمای جامع

سیستم غیر خطی گستره وسیعی از سیستم‌های موجود در طبیعت را در بر می‌گیرد. در اکثر موارد، سیستم را می‌توان توسط یک مدل خطی ساده‌تر تقریب زد. در این موارد اثر غیر خطی بودن (Nonlinearity) کوچک است و می‌توان از آن صرف نظر کرد. اما این ایده همیشه درست نیست. بسته به نوع غیر خطی بودن، حالت‌های سیستم می‌توانند رفتار همگرا (Convergent)، واگرا (Divergent) یا نوسانی (Cyclic) داشته باشند. در این مطلب قصد داریم تا سیستم‌های غیر خطی را معرفی کنیم و به بیان روش‌های آنالیز پایداری بپردازیم. همچنین این مفاهیم را با ذکر چند مثال ساده، واضح‌تر خواهیم ساخت.

997696

یک مثال ساده از یک خودروی زیرآبی (Underwater Vehicle) را در نظر بگیرید که معادله آن به صورت زیر است:

x¨=ux˙x˙\ddot{x} = u - \dot{x} |\dot{x}|

عبارت x˙x˙\dot{x} |\dot{x}| نشان‌دهنده اصطکاک ناشی از آب است.

در این مورد، حالت خاصی را که در آن ورودی کنترل u=0u = 0 است بررسی می‌کنیم. پس معادله سیستم برابر است با:

x¨=x˙x˙\ddot{x} = - \dot{x} | \dot{x} |

این سیستم زمانی که X˙>0\dot{X}>0، دارای مشتق منفی و زمانی که X˙<0\dot{X}<0، دارای مشتق مثبت است. بنابراین نرخ رشد حالت‌های سیستم همواره کاهشی است. در نتیجه می‌توان گفت که سیستم همواره پایدار است. این نوع غیر خطی بودن را اتلافی (Dissipative) یا پایدار‌کننده (Stabilizing) می‌نامند.

اما همیشه به این صورت نیست و ممکن است غیر خطی بودن منجر به ناپایداری حالت‌های سیستم شود. سیستم زیر را در نظر بگیرید:

x¨=ux2\ddot{x} = u - x^2

در این سیستم اگر حالت‌های سیستم با یک مقدار مثبت شروع شوند، گرادیان آن‌ها منفی می‌شود و حالت‌های سیستم به صفر کاهش می‌یابند. اما اگر مقدار اولیه xx منفی باشد، گرادیان عددی منفی می‌شود و حالت سیستم به -\infty می‌رود. در این حالت نوع غیر خطی بودن، غیراتلافی (Non-dissipative) یا ناپایدارکننده (Destabilizing) نامیده می‌شود.

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

نوسان‌ساز وندر پل

یک سیستم غیر خطی مرتبه دوم که بسیار مورد بررسی قرار می‌گیرد، نوسان‌ساز وندر پل (Van der Pol Oscillator) است. فرم این سیستم به صورت زیر است:

X¨+u(X)X˙+V(X)=0\ddot{X} + u(X) \dot{X} + V(X) = 0

که در آن u(X)u(X) و V(X)V(X) می‌توانند به عنوان چسبندگی (Viscosity) و سختی (Stiffness) در نظر گرفته شوند. حالت خاصی را در نظر بگیرید که:

X¨μ(1X2)X˙+X=0\ddot{X} - \mu (1- X^2) \dot{X} + X = 0

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

X2˙=μ(1X12)X2X1\dot{X_2} = \mu (1- X_1^2) X_2 - X_1

در معادله بالا μ\mu پارامتر میرایی (Damping) است. یک راه متداول برای به تصویر کشیدن رفتار یک سیستم غیر خطی، رسم یک حالت بر حسب حالت‌های دیگر است و معمولا یک متغیر بر حسب مشتق آن رسم می‌شود. به دیاگرام ترسیم‌شده پرتره فاز می‌گویند. قطعه کد زیر، نحوه رسم حالت‌های سیستم بالا در MATLAB را نشان می‌دهد.

1clc
2close all
3clear all
4
5mu = 1;
6
7sys_vanderPol = @(t,X)[X(2);mu*(1-X(1)^2)*X(2) - X(1)];
8
9X0 = [1;0];
10time_span = 0:.1:40;
11
12[time,states] = ode45(sys_vanderPol,time_span,X0);
1figure;
2subplot(2,1,1)
3plot(time,states(:,1))
4xlabel('time')
5ylabel('X')
6subplot(2,1,2)
7plot(time,states(:,2))
8xlabel('time')
9ylabel('dX/dt')

شکل زیر حالت‌های سیستم را نشان می‌دهد.

رسم حالت‌ها در یک معادله وندر پل
رسم حالت‌ها در یک معادله وندر پل

حال کد زیر پرتره فاز سیستم را رسم می‌کند.

1figure;
2plot(states(:,1),states(:,2));
3xlabel('X');
4ylabel('dX/dt');

پرتره فاز سیستم به صورت زیر است.

پرتره فاز سیستم وندر پل
پرتره فاز سیستم وندر پل

همان‌طور که از شکل بالا مشخص است، علی‌رغم شرایط اولیه سیستم، حالت‌ها به یک چرخه نوسانی همگرا می‌شوند. به این چرخه نوسانی، چرخه حدی (Limit Cycle) می‌گویند. بنابراین علاوه بر این‌که حالت‌ها می‌توانند به یک مقدار ثابت همگرا شوند و یا به بی‌نهایت واگرا شوند، این امکان نیز وجود دارد که حالت‌ها رفتاری دوره‌ای از خود نشان دهند که به آن چرخه حدی می‌گویند. طبقه خاصی از سیستم‌های غیر خطی که چنین رفتار نوسانی را از خود نشان می‌دهند، سیستم‌های لینارد (Lienard) نامیده می‌شوند.

سیستم‌های لینارد

معادله دیفرانسیل زیر را در نظر بگیرید:

y1=f1(t,X)y_1 = f_1(t,X)

y2=f2(t,X)y_2 = f_2(t,X)

اگر به ازای چند f1f_1 و f2f_2 داشته باشیم:

f1t+f2t0\frac{\partial f_1}{\partial t} + \frac{\partial f_2}{\partial t} \neq 0

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

X¨+F(X)X˙+G(X)=0,\ddot{X} + F'(X) \dot{X} + G'(X) = 0,

تئوری لینارد برای چرخه حدی

تئوری لینارد بیان می‌کند که اگر دو تابع نرم و فرد FF و GG' داشته باشیم، به نحوی که برای X>0X>0، G(X)>0G'(X)>0 و همچنین تابع FF دقیقا دارای سه صفر در 00 و aa و a-a باشد و F(0)<0F’(0) < 0 و F(X)0F’(X) ≥ 0 برای X>aX > a برقرار باشند و نیز F(X) F(X) → ∞ اگر X X →∞، آنگاه سیستم Lienard متناظر با آن دقیقا دارای یک چرخه حدی است و این چرخه پایدار است.

توجه کنید که عبارت بالا بیان می‌کند که مشتق تابع FF اگر XX بزرگتر از صفر باشد، افزایش می‌یابد و اگر XX بزرگتر از aa باشد، کاهش می‌یابد. این شروط منجر به سیستمی می‌شوند که در آن حالت‌های سیستم از صفر دور شده و به جای همگرایی به یک مقدار ثابت، به یک مسیر دورانی (Orbit) نوسانی همگرا می‌شوند. اگر FF و GG به صورت زیر انتخاب شوند:

F(X)=μ(X33X)F(X) = \mu \left( \frac{X^3}{3}-X \right)

G(X)=X22,G(X) = \frac{X^2}{2},

X¨+F(X)X˙+G(X)=0\ddot{X} + F'(X) \dot{X} + G'(X) = 0

آن‌گاه سیستم به صورت زیر نوشته خواهد شد:

X¨μ(1X2)X˙+G(X)=0\ddot{X} - \mu \left( 1 - X^2 \right) \dot{X} + G'(X) = 0

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

تکنیک‌های آنالیز پایداری سیستم‌ غیر خطی

سیستم غیر خطی در دو مد، رفتار پایداری از خود نشان می‌دهند. این سیستم‌ها ممکن است یا یک نقطه تعادل (Equilibirum) داشته باشند که حالت‌های سیستم به آن همگرا شوند و یا یک چرخه حدی داشته باشند که حالت‌ها رفتاری نوسانی را تجربه کنند.

نقطه تعادل

سیستمی را که توسط معادله زیر توصیف می‌شود در نظر بگیرید:

X˙=f(X)\dot{X} = f(X)

یک نقطه XeqX_{eq} نقطه تعادل نامیده می‌شود اگر شرط زیر برقرار باشد:

f(Xeq)=0f( X_{eq}) = 0

این نقطه تعادل باید پایدار باشد. اگر در همسایگی XeqX_{eq}، جهت مشتق از آن دور شود، حالت‌ها حرکت خواهند کرد.

چرخه حدی

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

r˙=(r21)\dot{r} = (r^2-1)

θ˙=1\dot{\theta} = 1

در این سیستم حالت‌های سیستم تا r=1r=1 همگرا می‌شوند، اما برای r>1r>1، حالت‌ها به بی‌نهایت می‌روند.

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

  1. حل معادلات دیفرانسیل به روش عددی (شبیه‌سازی) یا به روش تحلیلی و مطالعه رفتار جواب.
  2. خطی‌سازی سیستم غیر خطی در تعدادی از نقاط و مطالعه رفتار آن‌ها. این روش مستعد خطای خطی‌سازی است.
  3. تعریف یک تابع لیاپانوف و مطالعه رفتار خود تابع.

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

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

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

  1. ابتدا تابع LL را به نحوی تعریف می‌کنیم که در تمام نقاط فضا به غیر از نقطه تعادل مثبت معین باشد.
  2. مشتق این تابع را نسبت به زمان محاسبه می‌کنیم و کنترل می‌کنیم که تابع دارای کدام‌یک از مشخصه‌های زیر است:
    • اگر مشتق این تابع همواره منفی باشد، آن‌گاه حالت‌ها به صورت مجانبی (Asymptotically) به نقطه تعادل میل خواهند کرد.
    • اگر مشتق این تابع برای تمام X>aX>a به ازای یک aa معین منفی شود، آن‌گاه aa نشان‌دهنده مرز پایداری است. به عبارتی دیگر، بعد از زمان بی‌نهایت حالت‌های سیستم توسط aa برای همیشه محدود می‌شوند.

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

حال از روش لیاپانوف برای مطالعه رفتار یک سیستم گسسته در زمان مانند یک نوسان‌ساز وندر پل استفاده می‌کنیم. سیستم زیر را در نظر بگیرید:

X1˙=X2\dot{X_1} = X_2

X2˙=μ(1X12)X2X1\dot{X_2} = \mu (1- X_1^2) X_2 - X_1

تابع لیاپانوفی را به صورت زیر در نظر می گیریم:

L(X1,X2)=12X1TX1+12X2TX2L(X_1,X_2) = \frac{1}{2} X_1^T X_1 + \frac{1}{2} X_2^T X_2

حال از این معادله نسبت به زمان مشتق می‌گیریم:

Lt=X1TX˙1+X2TX˙2\frac{ \partial L }{ \partial t} = X_1^T \dot{X}_1 + X_2^T \dot{X}_2

با جایگذاری معادلات حالت در مشتق لیاپانوف داریم:

Lt=X1TX2+X2T(μ(1X12)X2X1)\frac{ \partial L }{ \partial t} = X_1^T X_2 + X_2^T \left( \mu (1- X_1^2) X_2 - X_1 \right)

Lt=μX2T(1X12)X2\frac{ \partial L }{ \partial t} = \mu X_2^T (1- X_1^2) X_2

توجه کنید که μ=0\mu=0 یک حالت خاص است که در آن تابع لیاپانوف تغییری نمی‌کند و سیستم یک حرکت تناوبی ثابت را از خود نشان می‌دهد. بنابراین بر روی دو حالت دیگر یعنی μ>0\mu >0 و μ<0\mu < 0 تمرکز می‌کنیم:

حالت μ>0\mu > 0

در این حالت به سادگی μ\mu را با قدر مطلق آن جایگزین می‌کنیم. زیرا باعث هیچ تغییری در رفتار سیستم نخواهد شد:

Lt=μX2T(1X12)X2\frac{ \partial L }{ \partial t} = |\mu| X_2^T (1- X_1^2) X_2

بنابراین اگر 1X12<01 - X_1^2 < 0، تابع لیاپانوف کاهش خواهد یافت. پس متغیرهای حالت همیشه محدود خواهند ماند و اگر (1X12)>0(1 - X_1^2) > 0، حالت‌های سیستم به سمت X12=1X_1^2 = 1 خواهند رفت.

حالت μ<0\mu < 0

در این حالت μ\mu را با منفی قدر مطلق آن جایگزین خواهیم کرد. زیرا منجر به تغییر رفتار سیستم نخواهد شد. تابع لیاپانوف به صورت زیر در خواهد آمد:

Lt=μX2T(1X12)X2\frac{ \partial L }{ \partial t} = - |\mu| X_2^T (1- X_1^2) X_2

بنابراین اگر (1X12)>0(1 - X_1^2) > 0 باشد، حالت‌های سیستم به صفر همگرا خواهند شد و اگر (1X12)<0(1 - X_1^2) < 0 حالت‌های سیستم به بی‌نهایت خواهند رفت.

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

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

کدهای زیر نوسان‌ساز وندر پل را به ازای μ=1\mu=1 رسم می‌کند:

1clc
2close all
3clear all
4
5X1_0 = -4:.45:4;
6X2_0 = -4:.45:4;
7
8[X1_0,X2_0] = meshgrid(X1_0,X2_0);
9mu = 1;
10sys_vanderPol = @(t,X)[X(2);mu*(1-X(1)^2)*X(2) - X(1)];
11
12X0 = [X1_0(1);X2_0(1)];
13time_span = 0:.1:40;
14
15[time,states] = ode45(sys_vanderPol,time_span,X0);
16[time,states1] = ode45(sys_vanderPol,time_span,[.2;0]);
17
18dX1  = X2_0;
19dX2  = (1-X1_0.^2).*X2_0 - X1_0;
20dX1 = dX1./sqrt(dX1.^2+dX2.^2);
21dX2 = dX2./sqrt(dX1.^2+dX2.^2);
1figure;
2subplot(1,2,1)
3plot(states(:,1),states(:,2),'linewidth',3)
4hold on
5quiver(X1_0,X2_0,dX1,dX2,.5,'r')
6axis([-4 4 -4 4])
7xlabel('X');
8ylabel('dX/dt');
9subplot(1,2,2)
10plot(states1(:,1),states1(:,2),'linewidth',3)
11hold on
12quiver(X1_0,X2_0,dX1,dX2,.5,'r')
13axis([-4 4 -4 4])
14xlabel('X');
15ylabel('dX/dt');

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

چرخه حدی پایدار در سیستم غیر خطی
چرخه حدی پایدار در سیستم غیر خطی

کد زیر نیز پرتره فاز را برای μ=1\mu=-1 رسم خواهد کرد.

1clc
2close all
3clear all
4
5X1_0 = -4:.45:4;
6X2_0 = -4:.45:4;
7
8[X1_0,X2_0] = meshgrid(X1_0,X2_0);
9mu = -1;
10sys_vanderPol = @(t,X)[X(2);mu*(1-X(1)^2)*X(2) - X(1)];
11
12time_span = 0:.1:10;
13
14[time,states] = ode45(sys_vanderPol,time_span,[2.1;0]);
15[time,states1] = ode45(sys_vanderPol,time_span,[2;0]);
16
17dX1  = X2_0;
18dX2  = (1-X1_0.^2).*X2_0 - X1_0;
19dX1 = dX1./sqrt(dX1.^2+dX2.^2);
20dX2 = dX2./sqrt(dX1.^2+dX2.^2);
1figure;
2subplot(1,2,1)
3plot(states(:,1),states(:,2),'linewidth',3)
4hold on
5quiver(X1_0,X2_0,dX1,dX2,.5,'r')
6axis([-4 4 -4 4])
7xlabel('X');
8ylabel('dX/dt');
9subplot(1,2,2)
10plot(states1(:,1),states1(:,2),'linewidth',3)
11hold on
12quiver(X1_0,X2_0,dX1,dX2,.5,'r')
13axis([-4 4 -4 4])
14xlabel('X');
15ylabel('dX/dt');

شکل زیر پرتره فاز در این حالت را نشان می‌دهد.

چرخه حدی ناپایدار در سیستم غیر خطی
چرخه حدی ناپایدار در سیستم غیر خطی

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

^^

بر اساس رای ۱۵ نفر
آیا این مطلب برای شما مفید بود؟
اگر بازخوردی درباره این مطلب دارید یا پرسشی دارید که بدون پاسخ مانده است، آن را از طریق بخش نظرات مطرح کنید.
منابع:
Vivek Yadav
۶ دیدگاه برای «سیستم غیر خطی — راهنمای جامع»

سلام . ممنون از آموزش کاربردی شما.
برای کدنویسی دینامیک غیرخطی در متلب پیشنهاد منبع میدید؟!

با سلام
وقتتون بخیر
برای خطی سازی سیستم های غیر خطی گسسته که متغیر با زمان هستند و پارامتر متغیر نیز دارند. فیلم یا مطلب اموزشی به من معرفی میکنید؟ ممنونم.

سلام بسیار مفید و قابل استفاده بود. در صورت امکان منبع مطالب بیان شده برای استفاده و مطالعه بیشتر را ذکر کنید.
بسیار ممنون

سلام خانم مهندس وقت بخیر
سوالی داشتم خدوتون : بنده دارم روی یک سیستم کار می کنم که در حالت پیوسته پایدار(قطب ها منفی) و در حالت گسسته نیز پایدار(قطب ها داخل دایره واحد) هستند و گسسته سازیمو به دو روش انجام دادم 1- اویلر 2 – رانگ کوتای مرتبه4 ولی متاسفانه سکل موج های خروجیم ناپایداره و سیگنال منترلیم ندسانات زیادی داره و دامنه آن نیز زیاد هستش تقریبا 500 یا شایدم بیشتر…. . به نظرتون برای پایداریش باید چیکار کنم؟؟
ممنونتون میشم راهنماییم کنید از طریق ایمیل .

با سلام
ممکن است گام زمانی شبیه‌سازی نزدیک به گام زمانی روش گسسته سازی انتخاب شده باشد. گام شبیه سازی را باید حداقل ۱۰ برابر کوچک‌تر از گام گسسته‌سازی انتخاب کرده و حتی مقدار ۱۰۰ برابر کوچک‌تر را نیز امتحان کنید.

سلام
خیلی مفید، مختصر و آموزنده بود.
تشکر

نظر شما چیست؟

نشانی ایمیل شما منتشر نخواهد شد. بخش‌های موردنیاز علامت‌گذاری شده‌اند *