شبیه سازی مدار در متلب – راهنمای کاربردی (+ دانلود فیلم آموزش رایگان)

در آموزش‌های پیشین مجله فرادرس، با برنامه‌نویسی در متلب و اصول شبیه‌سازی در سیمولینک آشنا شدیم. در این آموزش، شبیه سازی مدار در متلب و سیمولینک را در قالب مثال یک مدار RLC بیان می‌کنیم.

مدل‌سازی مدار RLC

یک مدار RLC در شکل ۱ نشان داده شده است.

ابتدا معادلات مدار RLC را به دست می‌آوریم. در مدار مورد نظر، جریان $$i$$ رابطه مستقیمی با افت ولتاژ‌ $$v _L$$ سلف دارد؛

فیلم آموزش فیزیک 2 دانشگاه – جامع و با نکات کاربردی در فرادرس

کلیک کنید

بدین صورت که مشتق ولتاژ $$v_C$$ متناسب با جریان $$i$$ گذرنده از آن است. بنابراین، می‌توانیم دو معادله دیفرانسیل مرتبه اول و سپس با اعمال قانون ولتاژ کیرشهف (KVL)، معادله ولتاژ حلقه را بنویسیم.

$$\large L \frac {d i } { d t } = v _ L \\ \large C \frac { d v _ C } { d t } = i \\ \large - v _ s + R i + v _L + v _ C = 0$$

دو متغیر $$i$$ و $$v _ C$$ متغیرهای حالت نامیده می‌شوند. معادلات بالا، به عنوان معادلات جبری دینامیکی (DAE) شناخته می‌شوند. با جایگذاری مقدار $$v _L$$ در معادله جبری، به دو معادله دیفرانسیل معمولی (ODE) خواهیم رسید:

$$\large L \frac { d i } { d t } = v _ s - R i - v _ C \\ \large C \frac { d v _c } { d t } = i$$

در معادلات بالا، مشتق متغیرهای حالت برحسب خود حالت‌ها و ورودی $$v _ s$$ بیان شده‌ است.

توصیف عمومی سیستم بالا به صورت زیر است:

$$\large \dot {x }  = f ( x , u )$$

که در آن، $$u$$ ورودی و $$x$$ بردار متغیرهای حالت است.

در مدار RLC بالا، $$x = [ i , v _ C ] ^ T$$ و $$u = v _ s$$. مدار RLC یک سیستم خطی است و می‌توانیم مدل فضای حالت خطی زیر را برای آن بنویسیم:

$$\large \dot { x } = \underbrace { \left [ \begin {array} {cc}{ \frac { - R } { L } } & { \frac { - 1 } { L } } \\ { \frac { 1 } { C } } & { 0 } \end {array} \right ] } _ { A } x + \underbrace { \left [ \begin {array} { c } { \frac { 1 } { L } } \\ { 0 } \end {array} \right ] } _ { B } \underbrace { v _ { s } } _ { u }$$

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

$$\large y = C x + D u$$

که در آن، $$C = I$$ یک ماتریس همانی است و $$D = [0,0]^ T$$.

اگر فقط جریان را به عنوان خروجی در نظر بگیریم، داریم:

$$\large y = [ 1 \;\;\;\; 0 ] x$$

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

برای مقایسه، شبیه‌سازی را در پایتون نیز انجام می‌دهیم. کد پایتون زیر برای شبیه‌سازی مدار RLC نوشته شده است. در این برنامه، ابتدا تابعی برای محاسبه مشتق‌های $$dx/dt$$ تعریف شده است.

فیلم آموزش برنامه نویسی پایتون Python – تکمیلی – بخش یکم در فرادرس

کلیک کنید

ورودی‌های این تابع پارامترهای RLC، بردار متغیر حالت $$x$$ و منبع ولتاژ $$v _ s$$ هستند. در ادامه برنامه، روش ذوزنقه‌ای رایج برای انتگرال‌گیری عددی نوشته شده است. در هر گام $$k$$، مشتق‌های $$dx/dt$$ با فراخوانی تابع محاسبه می‌شوند. بر اساس مشتق، حالت در گام‌های بعدی با استفاده از روش اویلر رو به جلو تخمین زده می‌شود. بر اساس این تخمین، مشتق در گام $$k + 1$$ محاسبه می‌گردد. با استفاده از دو مشتق، حالت $$x _ {k+1}$$ در گام $$k + 1$$ محاسبه خواهد شد.

1import math,pylab
2# define the dynamic equation to compute di/dt
3# R*i + L. di/dt +vc = v # c
4dvc/dt = i def fun_dxdt(R, L, C, vs, x):
5dx_dt = [-R/L*x[0]-x[1]/L+vs/L, x[0]/C]
6return dx_dt
7# Trapezoidal method to conduct numerical integration
8# step 1, initial condition
9# voltage is a dc voltage. for all time being, voltage is 1V.
10step_size = 0.001
11n_steps = 1000
12v = 1
13R = 0.1
14L = 0.01
15C =0.001
16x_data =[]
17v_data =[]
18x = [0, 0]
19x1 =[0, 0]
20x2 =[0, 0]
21for k in range(n_steps):
22v_data.append(v)
23x_data.append(x[:])
24x = x2
25# compute current
26dx_dt = fun_dxdt(R, L, C, v, x)
27for i in range(2):
28x1[i] = x[i] + dx_dt[i]*step_size
29dx_dt_est = fun_didt(R, L, C, v, x1)
30for i in range(2):
31x2[i] = x[i] + 0.5*(dx_dt[i]+dx_dt_est[i])*step_size
32tt = [k*step_size for k in range(n_steps)]
33pylab.plot(tt, x_data)
34pylab.xlabel('time (sec)');
35pylab.show()

شکل ۲ نتایج شبیه‌سازی برنامه پایتون بالا را نشان می‌دهد.

شبیه سازی مدار در متلب با حل کننده ODE

حل‌کننده‌های ODE متلب محیطی را برای انتگرال‌گیری عددی فراهم می‌کنند. تنها کافی است $$f ( x , u )$$ را تعریف کرده و حل‌کننده را فراخوانی کنیم.

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

کلیک کنید

در یک ام‌فایل (m-file) به نام fun_RLC.m، تابعی با نام fun_RLC را برای محاسبه $$f ( xx , u )$$ تعریف می‌کنیم. برای استفاده از حل‌کننده‌های ODE، ابتدا $$u$$ را تعریف می‌کنیم. در اینجا، $$u = 1$$ را به عنوان یک ورودی پله در نظر می‌گیریم تا پاسخ پله را مشاهده کنیم.

1function xdot = fun_RLC(t,x,par)
2xdot = [ -par.R/par.L, -1/par.L;1/par.C, 0]*x + [1/par.L; 0];

برنامه اصلی فراخوانی حل‌کننده ODE به صورت زیر است:

1par.R = 0.1; par.L = 0.01; par.C = 0.001;
2[t,y]=ode45(@(t,x) fun_RLC(t,x,par),[0 1],[0 0]);
3plot(t,y);

دقت کنید که در حل‌کننده‌های ODE باید نام تابع را برای محاسبه مشتق بردار حالت نوشت. ورودی‌های دوم و سوم تابع $$\text{ode45}$$ مدت زمان شبیه‌سازی و متغیرهای حالت اولیه هستند. همان‌طور که می‌بینیم، نتایج شبیه‌سازی مشابه شکل ۲ است.

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

در این بخش، دو روش را برای شبیه‌سازی مدار در سیمولینک بیان می‌کنیم. البته شبیه‌سازی را می‌توان با کمک کتابخانه‌هایی مانند Simscape Electrical در سیمولینک نیز انجام داد که در این آموزش به آن نمی‌پردازیم.

تشکیل مدل با تابع توکار Embedded Function

گام مهم فرایند تشکیل مدل شناسایی بردار متغیرهای حالت $$x$$ و محاسبه مشتق $$\dot {x } = f ( x , u )$$ از بردار ورودی $$u$$ و خودِ $$x$$ است. پس از آن، حالت و مشتق آن با یک انتگرال‌گیر با هم مرتبط می‌شوند. با تنظیم حالت اولیه $$x ( 0 )$$ در انتگرال‌گیر، مدل آماده اجرا خواهد بود.

شکل ۳ بلوک‌های تشکیل دهنده یک مدل دینامیکی را نشان می‌دهند که به صورت $$\dot { x } = f ( x , u )$$ توصیف می‌شود. از یک انتگرال‌گیر نیز برای برقراری ارتباط حالت و مشتق آن استفاده شده است. ورودی انتگرال‌گیر باید $$f ( x , u )$$ باشد، در حالی که خروجی آن، $$x$$ است. $$f (x , u )$$ با استفاده از یک تابع توکار (Embedded Function) متلب محاسبه می‌شود. ورودی تابع یک بردار شامل منبع ولتاژ $$v _ s$$ و متغیرهای حالت $$x$$ است. $$v _ s$$ و $$x$$ به بلوک Mux وارد می‌شوند و یک بردار را تولید می‌کنند. درون تابع fcn_RLC، $$v _ s$$ و $$x$$ از هم تفکیک شده‌اند. توجه کنید که ورودی دیگری به نام par نیز وجود دارد که یک ساختار (Structure) شامل مقادیر $$R$$، $$L$$ و $$C$$ است.

1function x_dot = fcn_RLC(u, par)
2%#codegen
3vs = u(1);
4x = u(2:end);
5A = [ -par.R/par.L, -1/par.L;
61/par.C, 0];
7B = [ 1/par.L; 0];
8x_dot = A*x+B*vs;
9return

همان‌طور که در شکل ۴ نشان داده شده است، متلب par را به عنوان یک ورودی در نظر می‌گیرد. برای آنکه par به عنوان یک پارامترِ مشخص در نظر گرفته شود، باید از گزینه Edit Data استفاده کرده و آن را به عنوان یک پارارمتر انتخاب کنیم.

یک مورد مهم‌تر، تنظیم صحیح مقادیر اولیه انتگرال‌گیر،‌ مخصوصاً‌ اندازه $$x ( 0 )$$ است. در این مثال، $$x ( 0 )$$ باید یک بردار شامل مقدار جریان اولیه و ولتاژ اولیه باشد. بنابراین، مقدار بردار اولیه را [0,0] یا (zeros(2,1 قرار می‌دهیم.

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

$$x ( 0 )$$ " width="423" height="516">

تشکیل مدل بر اساس S-function

در رویکرد S-function، کل سیستم در یک بلوک تعریف می‌شود که در شکل ۶ نشان داده شده است. یک S-function قالب مشخص و ثابتی برای تعریف ابعاد سیستم، حالت اولیه، مشتق حالت و... دارد. برای مدار RLC مورد نظر، از قالب سیستم پیوسته استفاده می‌کنیم.

S-function سه مورد دارد که باید به آن‌ها توجه کنیم: (الف) مقداردهی اولیه، (ب) محاسبه مشتق و (ج) تعریف خروجی. تنظیمات بردار حالت اولیه ورودی، خروجی و حالت در ابتدای کار تنظیم می‌شوند. این کار با تابع mdlInitializeSizes انجام می‌شود. کار دوم، محاسبه مشتق است که با کمک تابع mdlDerivatives انجام می‌گیرد. کار سوم، تعریف خروجی است که با تابع mdlOutputs آن را انجام می‌دهیم. برنامه S-function به صورت زیر است:

1function [sys,x0,str,ts,simStateCompliance] = fun_RLC_s(t,x,u,flag)
2switch flag,
3% Initialization %
4case 0,
5[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
6% Derivatives %
7case 1,
8sys=mdlDerivatives(t,x,u);
9% Outputs %
10case 3,
11sys=mdlOutputs(t,x,u);
12% other cases related to discrete systems
13case {2,4,9}
14sys=[];
15otherwise
16DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
17end
18function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
19sizes = simsizes;
20sizes.NumContStates = 2; % x is a vector of 2.
21sizes.NumDiscStates = 0;
22sizes.NumOutputs = 2;
23sizes.NumInputs = 1;
24sizes.DirFeedthrough = 0;
25sizes.NumSampleTimes = 1; % at least one sample time is needed
26sys = simsizes(sizes);
27% initialize the initial conditions
28x0 = [0;0];
29str = [];
30% initialize the array of sample times%
31ts = [0 0];
32simStateCompliance = 'UnknownSimState';
33% case 1: derivative
34function sys=mdlDerivatives(t,x,u)
35par.R = 0.1;
36par.L = 0.01;
37par.C = 0.001;
38A = [ -par.R/par.L, -1/par.L;
391/par.C, 0];
40B = [ 1/par.L; 0];
41x_dot = A*x+B*u;
42sys = x_dot;
43% case 3: output
44function sys=mdlOutputs(t,x,u)
45sys = x;

دستورهای متلب برای شبیه‌سازی سیستم‌های خطی

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

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

کلیک کنید

در ادامه، با در نظر گرفتن مدار RLC، مدل LTI آن را در قالب تابع تبدیل یا فضای حالت خواهیم ساخت.

تعریف مدل‌های خطی

اولین گام، تعریف مدل LTI مدار RLC است. در بخش‌های قبل، مدل این مدار را در فرم زیر به دست آوردیم:

$$\large \begin {align*} \dot { x } & = A x + B u \\ y & = C x + D u \end {align*}$$

با تعریف $$A$$، $$B$$، $$C$$ و $$D$$ این مدل فضای حالت را می‌توان، در یک خط، به صورت زیر تعریف کرد:

1sys1 = ss(A,B,C,D)

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

1sys2 = tf(sys1)

$$\text{sys2}$$ به فرم یک تابع تبدیل لاپلاس است.

یک رویکرد دیگر تعریف مدل در حوزه فرکانس است. بدین منظور، مدل امپدانس را برای هر عنصر مدار بررسی می‌کنیم. مدل امپدانس اصطلاحی است که نسبت نِمُو یا تغییر ولتاژ به نمو جریان $$\frac { \Delta V (s)} { \Delta I (s)}$$ را در حوزه فرکانس توصیف می‌کند. تکنیک مدل‌سازی امپدانس، یک تکنیک مدل‌سازی پرکاربرد در تحلیل مبدل‌های الکترونیک قدرت است. مدل امپدانس مقاومت همان $$R$$ و مدل امپدانس سلف و خازن نیز، به ترتیب، $$Ls$$ و $$\frac { 1 } { C s }$$ است.

بنابراین، می‌توانیم عبارت تغییر جریان را با قانون اهم به دست آوریم. تابع تبدیل $$\Delta V ( s)$$ به $$\Delta I ( s)$$ به صورت زیر است:

$$\large \frac { \Delta I (s)} { \Delta V _ s (s) } = \frac { 1 } { R + L s + 1/ (Cs)}$$

در یک سیستم خطی، هنگام توصیف تابع تبدیل، معمولاً از نوشتن $$\Delta$$ چشم‌پوشی می‌کنیم. همچنین $$s$$ داخل پرانتز را نمی‌نویسیم.

$$\large \frac { I } { V _ s } = \frac { 1 } { R + L s + 1/ (Cs)} \\ \large \frac { V _C } { V _ s } = \frac {I \times \frac {1 } {Cs} } {V _ s } \\ \large \frac { V _ L } { V _ s} = s L \frac { I } { V _ s}$$

توابع تبدیل بالا روابط ورودی/خروجی جریان، ولتاژ خازن و ولتاژ سلف در برابر ولتاژ‌ ورودی را تعریف می‌کنند. این توابع تبدیل را می‌توان با استفاده از عملگر لاپلاسِ s در متلب بیان کرد:

1s = tf('s');
2sys2 = 1/(R+ L*s + 1/(C*s));
3sys3 = sys2 *1/(C*s);
4sys4 = L*s*sys2;

خط دوم رابطه ورودی/خروجی ولتاژ و جریان منبع را تعریف می‌کند. خط سوم نیز رابطه ورودی/خروجی ولتاژ منبع و ولتاژ خازن را معرفی می‌کند. خط چهارم برای تعریف رابطه ورودی/خروجی بین ولتاژ منبع و ولتاژ‌ سلف است.

پاسخ زمانی

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

مثال ۱: پاسخ زمانی جریان $$i ( t)$$ را وقتی ولتاژ منبع یک پاسخ پله باشد.

حل: پاسخ پله را می‌توان با استفاده از دستور step متلب به دست آورد. برنامه متلب مورد نظر به صورت زیر است:

1R = 0.1; L = 0.01; C = 0.001;
2s = tf('s');
3sys2 = 1/(R+ L*s + 1/(C*s));
4step(sys2);

مثال ۲: پاسخ حوزه زمان جریان $$i(t)$$ را وقتی که ولتاژ‌ منبع از ۰ به یک ولتاژ‌سینوسی با فرکانس ۶ هرتز تغییر می‌کند، به دست آورید.

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

1T = 0:0.001:1.2;
2u = cos(37.7*T);
3lsim(sys2, u, T);

شکل ۷ نمودارهای مربوط بو دو دستور step و lsim را نشان می‌دهد.

تحلیل سیستم خطی

در یک سیستم تک‌ورودی-تک‌خروجی (SISO)، می‌توانیم تحلیل سیستم خطی را با استفاده از دستور bode متلب برای به دست آوردن نمودارهای بُد، دستور pole برای یافتن قطب‌ها، دستور zero برای یافتن صفرهای سیستم و pzmap برای پیدا کردن قطب‌ها و صفرها در فضای حقیقی-موهومی انجام دهیم.

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

1R1 = 0.1; R2 = 0.001;R3 = 1;
2L = 0.01;C = 0.001;
3s = tf('s');
4% current
5G1 = 1/(R1+ L*s + 1/(C*s));
6G2 = 1/(R2+ L*s + 1/(C*s));
7G3 = 1/(R3+ L*s + 1/(C*s));
8figure(1);
9bode(G1,G2,G3);
10grid;
11% capacitor voltage
12figure(2);
13h = bodeplot(G1/(C*s), G2/(C*s), G3/(C*s));
14grid;
15% Change units to Hz and make phase plot invisible
16setoptions(h,'FreqUnits','Hz','PhaseVisible','off');

$$I/V_s$$ و $$V_c/V_s$$" width="884" height="340">

نمودارهای بُد (اندازه) شکل ۸ مقدار پیک ۵۰۰ هرتز را نشان می‌دهند. این موضوع نشان دهنده وجود نوسان‌هایی با فرکانس ۵۰۰ هرتز است. اگر مقاومت کوچک‌تر باشد، مقدار این پیک محسوس‌تر خواهد بود. هرچه پیک نمودار اندازه بُد بیشتر باشد، نوسان در فرکانس ۵۰۰ هرتز میرایی کمتری خواهد داشت. با مقایسه شکل‌های ۷ و ۲ می‌توان دریافت که فرکانس‌ نوسان‌ها تقریباً ۵۰۰ هرتز است.

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

• مجموعه آموزش‌‌های برنامه‌‌ نویسی
• آموزش شبیه سازی با سیمیولینک Simulink
• مجموعه آموزش‌های برنامه‌نویسی متلب (MATLAB)
• گنجینه آموزش‌های شبیه سازی با سیمیولینک (Simulink)
• کاربرد تبدیل لاپلاس در تحلیل مدار — به زبان ساده
• رسم نمودار دو بعدی در متلب - از صفر تا صد
• دوقطبی در مدارهای الکتریکی — به زبان ساده

^^