پایداری سیگنال بزرگ سیستم قدرت — به زبان ساده

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

پایداری سیگنال بزرگ

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

فیلم آموزش بررسی سیستم های قدرت ۱ در فرادرس

کلیک کنید

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

در ادامه، معیار پایداری لیاپانوف را معرفی خواهیم کرد. سپس با استفاده از روش یا معیار مساحت مساوی (که مبتنی بر معیار پایداری لیاپانوف است)،‌ پایداری سیگنال بزرگ یا پایداری گذرای یک سیستم تک‌ماشین-شین بینهایت یا SMIB را بررسی می‌کنیم. نتایج تحلیل نیز با شبیه‌سازی حوزه زمان اعتبارسنجی شده است.

معیار پایداری لیاپانوف

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

$$\large \dot { x } = f ( x ) \; \; \; \; \; (1)$$

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

کلیک کنید

که در آن، $$x \in \mathbb { R }$$ و $$f : \mathbb { R } ^ n \to \mathbb { R } ^ n$$ است. قضیه پایداری لیاپانوف بیان می‌کند که این سیستم پایدار است، اگر و تنها اگر تابع $$V ( x ) \geq 0$$ که $$V : \mathbb { R } ^ n \to \mathbb { R }$$، وجود داشته باشد، به طوری که در آن، برای همه $$x$$های مسیر دینامیکی، مشتق $$V$$ کوچک‌تر یا مساوی با صفر شود ($$\dot{V} ( x ) \leq 0$$).

مثال ۱

با استفاده از قضیه پایداری لیاپانوف پایداری سیستم زیر را بررسی کنید.

$$\large \begin {align*} \dot {x } _ 1 & = - x _ 1 + g ( x _ 2 ) \\ \dot {x } _ 2 & = - x _ 2 + h ( x _ 1 ) \end {align*} \;\;\;\;\; (2)$$

که در آن:

$$\large \begin {align*} | g ( u ) | \leq | u | / 2 \\ | h ( u ) | \leq | u | / 2 \end {align*}$$

حل: تابع $$V ( x ) = \frac { 1 } { 2 } ( x _ 1 ^ 2 + x _ 2 ^ 2 )$$ را در نظر می‌گیریم که برای همه $$x$$ها $$V ( x ) \geq 0$$ است. اکنون، مشتق $$V ( x )$$ را محاسبه می‌کنیم:

$$\large \begin {align*} \dot { V } ( x ) & = x _ 1 \dot { x _ 1 } + x _ 2 \dot { x _ 2 } = x _ 1 ( - x _ 1 + g ( x _ 2 )) + x _ 2 ( - x _ 2 + h ( x _1)) \\ & = - x _ 1 ^ 2 - x _ 2 ^ 2 + x _ 1 g ( x _ 2 ) + x _ 2 h ( x _ 1 ) \\ & \leq - x _ 1 ^ 2 - x _ 2 ^2 + |x _ 1 | | g ( x _ 2 ) | + |x _ 2 | | h ( x _ 1 ) | \\ & \leq - x _ 1 ^ 2 - x _ 2 ^ 2 + | x _ 1 | \frac { | x _ 2 | } { 2 } + | x _ 2 | \frac {| x _ 1|} { 2 } \\ & = - x _ 1 ^ 2 - x _ 2 ^ 2 + | x _ 1 x _ 2 | \\ & \leq - \frac { 1 } { 2 } (x _ 1 ^ 2 + x _ 2 ^ 2 ) = - V ( x ) \;\;\;\;\; (3) \end {align*}$$

بنابراین، برای هر $$x$$ خواهیم داشت: $$\dot {V} ( x ) \leq 0$$. در نتیجه، سیستم بالا یک سیستم پایدار است.

دینامیک یک سیستم خطی تغییر ناپذیر با زمان (LTI) را می‌توان به صورت زیر نوشت:

$$\dot { x } = A x \;\;\;\;\; (4)$$

که در آن، $$A \in \mathbb {R} ^ { n \times n}$$ و $$x \in \mathbb{R}^n$$.

تابع $$V (x) = x ^ T P x$$ را در نظر بگیرید که در آن، $$P \in \mathbb {R} ^ { n \times n}$$ یک ماتریس نیمه‌معین مثبت (PSD) است؛ یعنی $$P$$ متقارن ($$P=P^T$$) بوده و هر مقدار ویژه آن بزرگ‌تر یا مساوی صفر است. این نوع ماتریس را با $$P \succeq 0$$ نشان می‌دهیم. در نتیجه، برای هر $$x$$، داریم: $$V (x) \geq 0$$.

گفته بالا را می‌توان با استفاده از تجزیه مقادیر ویژه اثبات کرد. ماتریس $$P$$ را می‌توانیم به صورت $$P = V ^ { - 1 } \Lambda V$$ تجزیه کنیم که در آن، $$\Lambda$$ یک ماتریس قطری با درایه‌های $$\lambda _ i \geq 0$$ $$i = 1 , \cdots , n )$$) و $$V$$ ماتریس بردار ویژه است. ماتریس $$V$$ متعامد است؛ یعنی $$V ^ T =V ^ {-1}$$.

$$\large V ( x ) = x ^ T V ^ T \Lambda V x = \tilde {x } ^ T \Lambda \tilde {x} = \sum _{ i = 1 }^{ n } { \lambda _ i \tilde {x } _ i ^ 2 \geq 0 \;\;\;\;\; ( 5 ) }$$

که در آن، $$\tilde {x} = V x$$ است. مشتق تابع لیاپانوف به صورت زیر محاسبه می‌شود:

$$\large \dot { V } (x ) = x ^ T P \dot { x } + \dot { x} ^ T P x = x ^T ( P A + A ^ T P ) x \;\;\;\;\; ( 6 )$$

اگر سیستم خطی پایدار باشد، آنگاه $$\dot {V} (x) \leq 0$$ است. این معادل با $$PA+A^T P \preceq 0$$ خواهد بود.

پایداری یا ناپایداری

قضیه پایداری لیاپانوف بیان می‌کند که سیستم پایدار است اگر $$V(x) \geq 0$$ بوده و برای همه $$x$$های مسیر $$\dot{V}(x) \leq 0$$ باشد. سیستم ناپایدار است، اگر برای همه $$x$$های روی مسیر، $$V(x) \lt0$$ بوده و $$\dot{V}(x) \leq 0$$ باشد.

روش مساحت مساوی

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

$$\large \begin {align*} \dot { \delta } & = \omega _ 0 ( \omega - 1 ) \\ \dot {\omega } & = \frac { 1 } { 2 H } ( P _ m - P _ e - D _ 1 ( \omega - 1 )) \end {align*} \;\;\;\;\; (7)$$

فیلم آموزش کاربرد فازی در سیستم های قدرت در فرادرس

کلیک کنید

تابع $$V ( \delta , \omega )$$ را تشکیل می‌دهیم. برای پایدار بودن، مشتق زمانی این تابع به ازای همه $$\delta$$ها و $$\omega$$ها باید کوچک‌تر یا مساوی صفر باشد. مشتق زیر را در نظر بگیرید:

$$\large \dot {V} ( \delta , \omega ) = - D_ 1 \omega _0 ^ 2 ( \omega - 1 ) ^ 2 = - D_ 1 \dot{\delta} ^ 2 \;\;\;\;\; (8 )$$

که به ازای هر $$(\delta , \omega)$$، نامساوی $$\dot{V} ( \delta , \omega ) \leq 0$$ برقرار است.

برای اظهار نظر درباره پایدار بودن یا نبودن سیستم، لازم است بررسی کنیم که برای همه $$(\delta , \omega)$$ها $$V (\delta , \omega ) \geq 0$$ باشد یا برای هر $$(\delta , \omega)$$ در مسیر، رابطه $$V ( \delta , \omega ) < 0$$ برقرار باشد.

$$\large \begin {aligned} V ( \delta, \omega ) & = \int _ { t _ { 0 } } ^ { t } \dot { V } ( \delta , \omega ) d t = \int _ { t _ { 0 } } ^ { t } - D _ {1 } ( \dot { \delta } ) ^ { 2 } d t \\ & = -\int _ { \delta _ { 0 } } ^ { \delta } D _ { 1 } \dot { \delta } d \delta \quad \\ & = -\int _ { \delta _ { 0 } } ^ { \delta } \left [ \omega _ { 0 } \left ( P _ { m } - P _ { e } \right ) - 2 H \omega _ { 0 } \dot { \omega } \right ] d \delta \\ & = -\omega _ { 0 } \int _ { \delta _ { 0 } } ^ { \delta } \left ( P _ { m } - P _ { e } \right ) d \delta + \int _ { \delta _ { 0 } } ^ { \delta } 2 H \frac { d \dot { \delta } } { d t } d \delta \\ & = -\omega _ { 0 } \int _ { \delta _ { 0 } } ^ { \delta } \left ( P _ { m } - P _ { e } \right ) d \delta + \int _ { \delta _ { 0 } } ^ { \delta } 2 H \dot { \delta } d \dot { \delta } \\ & = \omega _ { 0 } \int _ { \delta _ { 0 } } ^ { \delta } \left ( P _ {e } - P _ { m } \right ) d \delta + \left . H \left ( \dot { \delta } ^ { 2 } \right ) \right | _ { \delta _ { 0 } } ^ { \delta } \end {aligned} \;\;\;\;\; (9)$$

در شرایط اولیه $$\delta = \delta _ 0$$، رابطه $$\dot {\delta} = 0$$ را داریم. در یک نقطه مشخص از مسیر، وقتی سرعت به ۱ پریونیت می‌رسد (یا $$\dot{\delta} =0$$)، تنها لازم است $$\int _ {\delta _ 0 } ^ { \delta} { ( P _ e - P _ m )} d\delta < 0$$ را بررسی کنیم. اگر چنین باشد، آنگاه سیستم ناپایدار است. این نقطه از مسیر، متناظر با یک زاویه ماکزیمم $$\delta _ {\text{max}}$$ یا یک زاویه مینیمم $$\delta _ \text{min}$$ است.

بنابراین، اگر زاویه $$\delta$$ روتور افزایشی باشد، با معیار مساحت مساوی عبارت زیر محاسبه می‌شود:

$$\large \int _ { \delta _ 0 } ^ { \delta _ {\mathrm{max}}} ( P _ e - P_ m ) d \delta$$

مثال ۲

در یک سیستم SMIB، شین ترمینال ژنراتور در لحظه $$t_0$$ دچار خطای سه فاز به زمین شده است. زاویه اولیه روتور $$\delta _ 0$$ است. خطا در زمان $$t _ 1 = t_ 0 + \Delta  t$$ که زاویه روتور $$\delta _ 1$$ است، رفع می‌شود. پایداری سیستم را بررسی کنید. فرض کنید توان مکانیکی ثابت نگه داشته می‌شود. اگر سیستم پایدار است، مقدار $$\delta _ \text{max}$$ را تعیین کنید.

حل: $$\int _ { \delta _ 0 } ^ { \delta _ {\mathrm{max}}} ( P _ e - P_ m ) d \delta$$ را بررسی خواهیم کرد. این عبارت را می‌توان به دو بخش تقسیم کرد:

$$\large \int _ { \delta _ 0 } ^ { \delta _ {\mathrm{max}}} ( P _ e - P_ m ) d \delta = - \underbrace{ \int _ { \delta _ 0 } ^ { \delta _ 1 } ( P _ m - P_ e ) d \delta }_{A_1} + \underbrace{ \int _ { \delta _ 1 } ^ { \delta _ \text{max }} ( P _ e - P_ m ) d \delta }_{ A _ 2}. \;\;\;\;\; (10)$$

از $$t_0$$ تا $$t_1$$ خطای اتصال کوتاه رخ می‌دهد و ولتاژ ترمینال ژنراتور 0 و در نتیجه، توان حقیقی آن $$P_e = 0$$ است. بنابراین، از $$t_0$$ تا $$t_1$$، توان خالصی که در ژنراتور وجود دارد، سبب شتاب گرفتن آن می‌شود. $$A_1$$ مساحت شتاب نامیده می‌شود. $$A_2$$ نیز مساحت کاهش سرعت نام دارد، زیرا بعد از $$t_1$$، توان الکتریکی $$P_e$$ بازیابی می‌شود و $$P _e > P_m$$ خواهد شد.

$$\large A_1 = P_m ( \delta _1 - \delta _ 0 )$$

حداکثر مقدار $$\delta _ \text{max}$$، می‌تواند $$\pi - \delta _ 0$$ باشد. از آنجایی که اگر $$\delta > \delta _\text{max}$$ باشد، $$P_e< P_m$$ خواهد بود، ژنراتور شتاب می‌گیرد. مساحت شتاب بزرگ‌تر از $$A_1$$ است. بنابراین، حداکثر $$A_2$$ به صورت زیر خواهد بود:

$$\large \int _ { \delta _ 1 } ^ { \pi - \delta _ 0 } \left ( \frac {E V _ \infty } { X } \sin (\delta) - P_ m \right ) d \delta$$

اگر $$A_2 \geq A_ 1$$ باشد، سیستم پایدار است. در غیر این صورت، سیستم پایدار نیست. شکل ۱ این دو مساحت را نشان می‌دهد.

اگر سیستم پایدار، و $$\delta _ 0$$ و $$\delta _ 1$$ مشخص باشند، می‌توانیم $$\delta _ \text{max}$$ را با استفاده از معادله زیر بیابیم:

$$\large A_ 1 = P_ m ( \delta _ 1 - \delta _ 0 ) = A_ 2 = \int _ { \delta _ 1 } ^ {\delta _ \text{max} } \left ( \frac {E V _ \infty } { X } \sin (\delta) - P_ m \right ) d \delta \;\;\;\;\; (11)$$

مثال ۳

در یک سیستم SMIB، شین ترمینال ژنراتور در لحظه $$t_0$$ دچار خطای سه فاز به زمین شده است. در این لحظه، زاویه روتور $$\delta _ 0$$ است. خطا در $$t_ 1 = t _ 0 + \Delta t$$ برطرف می‌شود. حداکثر مقدار $$\Delta t$$ (یا زمان رفع بحرانی) چقدر باشد تا سیستم پایدار گردد؟

حل: زمان رفع بحرانی را می‌توان ابتدا از $$A_ 1 = A_2$$ و یافتن $$\delta _ 1$$ به دست آورد. علاوه بر این، $$\delta _ 1$$ را برحسب $$\Delta t$$ می‌نویسیم.

در طول $$t_ 0$$ تا $$t_ 1$$، $$P_ e = 0$$ است. بنابراین، اگر از اصطکاک مکانیکی صرف‌نظر کنیم ($$D_ 1 = 0$$)، معادله نوسان به صورت زیر در خواهد آمد:

$$\large \ddot { \delta } = \frac { \omega _ 0 } { 2 H } P _ m$$

توصیف حوزه زمان $$\delta _ 1$$ به صورت زیر است:

$$\large \delta _ 1 = \delta _ 0 + \frac { \omega_ 0 } { 2 H } P _ m (\Delta t ) ^ 2$$

در سیستمی با $$H = 5 \; \text{s}$$، $$D_ 1 = 0$$، $$P_m=1$$ و $$P_e = 2 \sin \delta$$، زاویه اولیه $$\delta _0$$ روتور (0٫5236 رادیان) را می‌توان با توجه به نقطه کار به عنوان یک نقطه تعادل محاسبه کرد.

$$\large \begin {align*} 0 & = \dot { \delta } = \omega _0 ( \omega - 1 ) \implies \omega = 1 \\ 0 & = \dot {\omega} = \frac { 1} { 2 H} (P_m - P_e -D_1 (\omega - 1)) =\frac{ 1 } {2 H} ( 1 - 2 \sin \delta ) \end {align*} \;\;\;\;\; (12)$$

زاویه بحرانی $$\delta _ 1 = 1.3886\; \text{rad}$$ با حل (۱۱) به دست می‌آید و زمان رفع بحرانی متناظر $$\Delta t = 0.2142\; \text{s}$$ است.

مثال ۴

اگر سیستم پایدار باشد، کمینه زاویه روتور $$\delta _ {\text{min}}$$ را به دست آورید.

حل: حداقل زاویه را نیز می‌توان با استفاده از روش مساحت مساوی یافت. نقطه شروع، همان زاویه روتور حداکثر است. با حل معادله زیر، مقدار $$\delta _ \text{min}$$ به دست خواهد آمد:

$$\large 0 = \int _{\delta _ \text{min} } ^ { \delta _ \text{max} } \left ( \frac { EV_ \infty } {X} \sin (\delta ) - P_ m \right ) d \delta$$

شکل ۲ دو ناحیه را نشان می‌دهد.

شبیه‌سازی در حوزه زمان

شبیه‌سازی زمانی منجر به اعتبارسنجی تحلیل نتایج از معیار ناحیه‌های برابر می‌شود.

فیلم آموزش دیگسایلنت – آنالیز و شبیه سازی سیستم های قدرت با DIgSILENT در فرادرس

کلیک کنید

برنامه متلب زیر، داده‌های زمانی مربوط به بردار متغیرهای حالت $$x = [\delta, \omega]^T$$ معادله نوسان و توان $$P_e$$ خروجی ژنراتور را تولید می‌کند. توجه کنید که روش انتگرال‌گیری عددی مورد استفاده، روش ذوزنقه‌ای است.

ورودی تابع، زمان رفع بحرانی $$tcr$$، گام زمانی $$h$$ و زمان پایان $$T\_ end$$ است.

1function [x, Pe]=SMIB_sim(tcr,h, T_end)
2omega_0= 377;
3parameter.H =5;
4parameter.D1 =0;
5delta0 = pi/6;
6E = 1; V=1; X=0.5;
7parameter.E = E;
8parameter.V = V;
9parameter.X = X;
10parameter.Pm = E*V/X*sin(delta0);
11parameter.tcr =tcr;
12x(:,1) = [delta0; 1];
13i=1;
14for t=0:h:T_end
15i= i+1;
16[dotx, Pe(i)] = SMIB(x(:,i-1), t, parameter);
17x1 = x(:,i-1) +dotx*h;
18[dotx1, Pe(i)] = SMIB(x1, t, parameter);
19x(:,i) = x(:,i-1) +(dotx+dotx1)/2*h;
20end
21return

این تابع، یک تابع SMIB دیگر را برای محاسبه $$\dot {x}$$ و $$P_e$$ در هر گام فراخوانی می‌کند. کد SMIB به صورت زیر است:

1function [x_dot, Pe] = SMIB(x, t, parameter)
2t0 = 1;
3x_dot = zeros(2,1);
4delta = x(1); omega = x(2);
5omega_0= 377;
6H = parameter.H;
7D1 = parameter.D1;
8Pm = parameter.Pm;
9tcr = parameter.tcr;
10E = parameter.E;
11V = parameter.V;
12X = parameter.X;
13if (t>t0 && t<t0+tcr)
14Pe = 0;
15else
16Pe = E*V/X*sin(delta);
17end
18x_dot(1) = omega_0*(omega-1);
19x_dot(2) = 1/(2*H)*(Pm - Pe -D1*(omega-1));

نتایج شبیه‌سازی زمانی یک سیستم SMIB در هنگام خطای سه فاز به زمین در شین ژنراتور، در شکل‌های ۳ تا ۵ نشان داده شده است.

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

1clear; clc;
2T_end = 3; h = 0.0002; T=0:h:T_end; n = length(T);
3[x, Pe]= SMIB_sim(0.2140,h, T_end); x2 = x; Pe2 = Pe;
4[x, Pe]= SMIB_sim(0.2142,h, T_end); x3 = x; Pe3 = Pe;
5[x, Pe]= SMIB_sim(0.2145,h, T_end); x4 = x; Pe4 = Pe;
6figure
7plot(T, [x2(1,1:n);x3(1,1:n); x4(1,1:n)],'LineWidth',2);
8ylabel('\delta'); grid on; xlabel('Time (s)'); xlim([0.9,2.2]);
9figure;
10plot(T, [x2(2,1:n);x3(2,1:n); x4(2,1:n)],'LineWidth',2);
11ylabel('\omega'); grid on; xlabel('Time (s)');xlim([0.9,2.2]);
12figure;
13plot(T, [Pe2(1:n);Pe3(1:n);Pe4(1:n)],'LineWidth',2);
14ylabel('P_eomega'); grid on; xlabel('Time (s)'); xlim([0.9,2.2]);

نتایج شبیه‌سازی تحلیل مساحت مساوی را تأیید می‌کنند. می‌توان دریافت که زمان رفع بحرانی ۰٫۲۱۴۲ ثانیه است.

$$\delta$$ (رادیان)‌" width="443" height="350">

$$\omega$$ (پریونیت)‌" width="445" height="350">

وقتی زمان رفع ۰٫۲۱۴۲ ثانیه باشد، از یک قاب زمانی طویل‌تر استفاده می‌شود. شکل ۶ نتایج شبیه‌سازی را برای ۲ ثانیه پس از خطا نشان می‌دهد.

توجه کنید که در زمان رفع بحرانی، بر اساس روش مساحت‌های برابر، $$\delta _ \text{max}$$ باید $$\pi - \delta _ 0 = \frac {5}{6} \pi = 2.618 \; \text{rad}$$ باشد. نتایج شبیه‌سازی حداکثر زاویه ۲٫۵۷۱۱ رادیان را نشان می‌دهد. خطا به دلیل خطای انتگرال‌گیری عددی است. اگر گام زمانی را کاهش دهیم، این خطا کاهش می‌یابد. برای مثال، اگر گام زمانی به $$10 ^ { -5}$$ ثانیه کاهش یابد، آنگاه $$\delta _ \text{max} = 2.595 \; \text{rad}$$ است.

$$P_e$$ (پریونیت)‌" width="462" height="350">

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

• مجموعه آموزش‌‌های مهندسی قدرت
• آموزش آنالیز پایداری و کنترل سیستم های قدرت با جعبه ابزارهای متلب
• مجموعه آموزش‌های مهندسی برق
• آموزش دینامیک سیستم های قدرت ۱
• مدل خط سه فاز — به زبان ساده
• کو انرژی چیست؟ — از صفر تا صد
• پروفایل بار چیست؟ — به زبان ساده

^^