پیوندهای مفید
اخبار و اعلانات
آمار
| تعداد نشریات | 43 |
| تعداد شمارهها | 1,791 |
| تعداد مقالات | 14,617 |
| تعداد مشاهده مقاله | 38,787,393 |
| تعداد دریافت فایل اصل مقاله | 15,100,387 |
طراحی ردیاب LPV در رباتهای سیار غیرهولونومیک در حضور خطای محرک و اشباع ورودی | ||
| هوش محاسباتی در مهندسی برق | ||
| مقاله 9، دوره 15، شماره 4، اسفند 1403، صفحه 115-130 اصل مقاله (915.42 K) | ||
| نوع مقاله: مقاله پژوهشی فارسی | ||
| شناسه دیجیتال (DOI): 10.22108/isee.2025.144574.1727 | ||
| نویسندگان | ||
| محمدکاظم شبانی گلویک1؛ محمد حسن آسمانی* 2؛ امیرحسین حامی3 | ||
| 1دانش آموخته کارشناسی ارشد مهندسی برق کنترل، دانشکدۀ مهندسی برق و کامپیوتر، دانشگاه شیراز، شیراز، ایران | ||
| 2استاد، دانشکدۀ مهندسی برق و کامپیوتر، دانشگاه شیراز، شیراز، ایران | ||
| 3دانشجوی کارشناسی ارشد رشتۀ مهندسی برق کنترل، دانشکدۀ مهندسی برق و کامپیوتر، دانشگاه شیراز، شیراز، ایران | ||
| چکیده | ||
| در سالهای اخیر، رباتهای چرخدار متحرک در بسیاری از زمینهها مانند پزشکی، کشاورزی و نظامی، توانستهاند جایگاهی ویژه پیدا کنند. این رباتها به دلیل طراحی ساده و هزینۀ کم تولید، به عنوان ابزارهایی کارآمد و مقرونبهصرفه برای انجام وظایف گوناگون در محیطهای مختلف شناخته میشوند. در این مقاله، یک راهبرد کنترلی جدید برای رباتهای سیار غیرهولونومیک مبتنی بر مدلهای LPV طراحی میشود. هدف از طراحی کنترلکنندۀ پیشنهادی تضمین پایداری دینامیکی سیستم و همچنین افزایش تابآوری آن در برابر خطاهای محرک و اشباع ورودی است که به طور معمول در سناریوهای عملیاتی رخ میدهند. فرایند طراحی با استفاده از حل یک دسته نامساویهای ماتریسی خطی (LMI) برای تعیین پارامترهای بهینۀ کنترلکننده انجام میشود. علاوه بر این، برای مقابله با خطاهای محرک، از فیلتر کالمن توسعهیافته (EKF) برای برآورد دقیق وضعیت سیستم و شناسایی تغییرات پارامترهای آن استفاده شده است. نتایج شبیهسازیها نشان میدهد کنترلکنندۀ پیشنهادی قادر به ردیابی دقیق مسیر مرجع است و در عین حال، میتواند پایداری و دقت خود را حتی در شرایطی که با خطاهای محرک و اشباع سیگنال کنترلی مواجه است، حفظ کند. | ||
| کلیدواژهها | ||
| ربات متحرک؛ ردیابی مسیر؛ کنترل تحملپذیر خطا؛ کنترلکنندۀ LPV؛ فیلتر کالمن توسعهیافته | ||
| اصل مقاله | ||
|
مقدمه
رباتهای سیار چرخدار، به دلیل قابلیت حرکت مستقل و انعطافپذیری زیاد، در صنایعی گوناگون مانند پزشکی، کشاورزی و حملونقل استفاده میشوند. یکی از چالشهای اصلی در کنترل رباتهای سیار تضمین پایداری و دقت مسیر حرکتی آنها در شرایط واقعی است. در بسیاری از کاربردها، این رباتها در محیطهای نامنظم یا پیچیده حرکت میکنند که ممکن است باعث بروز خطا در سیستمهای کنترلی شود ]1، 2[. در این راستا، طراحی کنترلکنندههای پیشرفتهای که بتوانند در مواجهه با این خطاها عملکردی مطلوب داشته باشند، ضروری به نظر میرسد ]3[.
روش LPV به عنوان یکی از روشهای پیشرفته و کارآمد در بهینهسازی عملکرد کنترل رباتهای سیار شناخته میشود؛ به طوری که این سیستمهای کنترلی قادر به تطبیق با تغییرات پارامترهای محیطی هستند و میتوانند مسیر خود را با دقت بیشتری دنبال کنند. در این روش، مدل سیستم به طور دینامیک و بر مبنای پارامترهای متغیر، بهروزرسانی میشود که در نتیجه، کنترلکننده را از انعطافپذیری بسیار زیادی برخوردار میکند ]4[.
یکی از مسائل مهم در کنترل مسیر رباتها مشکل اشباع ورودی در عملگرهای آنهاست ]5، 6[. معمولاً محرکهای رباتها محدودیتهایی از نظر توان عملیاتی دارند که این محدودیتها میتوانند دقت اجرای دستورات کنترلی را تحت تأثیر قرار دهند. برای مقابله با این چالش، میتوان اشباع ورودی را به عنوان یک محدودیت در طراحی کنترلکنندههای LPV در نظر گرفت ]7، 8[ و از روشهای بهینهسازی برای تنظیم دقیق پارامترهای کنترلکننده استفاده کرد. این رویکرد موجب جلوگیری از ناپایداری سیستم در شرایط مختلف عملیاتی خواهد شد ]8، 9[.
در بسیاری از پژوهشهای پیشین، طراحی کنترلکنندههای LPV برای رباتهای سیار، بدون لحاظکردن خطاهای ناشی از نقایص در عملکرد محرکها، انجام شده است. علاوه بر این، در برخی از این پژوهشها، فرض بر این بوده است که عملگرهای ربات قادر به تولید و اعمال ورودیهای کنترلی بدون هیچگونه محدودیت یا موانعی هستند ]7، 8[، در حالی که محدودیتهای سختافزاری عملگرها باعث اشباع ورودی و تأثیر آن بر عملکرد کنترلکننده میشوند. با توجه به نکات بیانشده، یک کنترلکنندۀ LPV نوآورانه ارائه خواهد شد که علاوه بر تضمین پایداری و ردیابی مسیر، قابلیت تحمل خطاهای محرکها را در حضور اشباع عملگرها، با استفاده از فیلتر کالمن توسعهیافته (EKF ) دارد ]10، 11[.
در این مقاله، یک راهبرد کنترلی جدید مبتنی بر مدلهای LPV برای رباتهای سیار غیرهولونومیک طراحی شده و توسعه یافته است که بتواند علاوه بر تضمین پایداری دینامیکی مسیر، عملکردی مقاوم در برابر انواع خطاهای ناشی از اختلالات در محرکها فراهم آورد. در این راستا، ابتدا یک روش جدید برای مدلسازی دقیق ربات سیار به صورت LPV ارائه میشود. سپس، شرایط پایداری دینامیک خطای ردیابی در حضور قیود ناحیۀ پایداری LMI ارائه میشود. در ادامه، از فیلتر کالمن توسعهیافته برای جبران خطاهای محرک استفاده میشود. علاوه بر این، اثرات منفی اشباع عملگر ورودی بر دقت و پایداری عملکرد کنترلکننده بررسی و راهکارهای مقابله با این اثرات ارائه شده است. در واقع، جنبههای نوآورانۀ این پژوهش شامل ارائۀ روشی جدید برای مدلسازی LPV رباتهای سیار، ارزیابی عملکرد کنترلکنندۀ ردیاب LPV در شرایط مختلف مانند خطای محرک، اشباع ورودی و تغییرات دینامیکی پارامترها هستند. با استفاده از نتایج شبیهسازی، نشان داده خواهد شد کنترلکنندۀ LPV طراحیشده با دقت زیاد قادر به ردیابی مسیر مرجع و واکنش مؤثر به خطاهاست و در شرایط چالشبرانگیز محیط واقعی عملکردی مطلوب دارد.
این مقاله به صورت زیر ساختاربندی شده است. در بخش 2، تعاریف اولیه و مقدمات ریاضی از موضوع، مرور و بررسی خواهند شد. سپس در بخش 3، کنترلکنندۀ LPV برای ردیابی مسیر رباتهای سیار دارای خطای محرک و اشباع ورودی طراحی خواهد شد. پس از آن، در بخش 4، نیز نتایج شیبهسازی مقاله مشاهده خواهد شد. در نهایت، در بخش 5، نتیجهگیری ارائه خواهد شد.
تعاریف اولیه و مقدمات ریاضی
1-2- سیستمهای خطی تغییرپذیر با پارامتر
روش مدلسازی به صورت LPV یکی از رویکردهای پیشرفته در کنترل سیستمهای رباتیک است که مزایای زیادی برای رباتهای متحرک فراهم میآورد ]9[. در این روش، مدل سیستم به گونهای طراحی میشود که بتواند تغییرات پارامترها را در زمان واقعی دنبال کند و به طور پویا به آنها پاسخ دهد ]11[. یکی از مزایای این روش این است که سیستمهای LPV قابلیت مدلسازی و کنترل سیستمهای غیرخطی را بهراحتی فراهم میآورند ]12[. سیستمهایLPV قادر هستند به طور مؤثر تغییرات غیرخطی را در مدلهای سیستم مدیریت کنند، در حالی که از ساختار خطی برای طراحی کنترلکنندهها بهره میبرند. این ویژگی موجب کاهش پیچیدگیهای سیستم و در نتیجه بهبود عملکرد آن میشود.
برای مدلسازی سیستم غیرخطی به صورت LPV میتوان از روش خطی سازی ژاکوبین استفاده کرد ]13[. یکی دیگر از روشهای مدلسازی تکنیک تبدیل حالت است که در مرجع ]14[بررسی شده است؛ به این صورت که ورودی کنترلی باعث حذف عبارتهای غیرخطی غیروابسته به پارامتر میشود. دیگر روشها، روش جایگزینی تابع ، روشهای ترکیبی با بهرهبردن از Linear Fits ]15[، روشهای تصادفی و روشهای خطیسازی بر اساس سرعت ]16[ هستند.
پایداری و کنترل سیستمهای LPV
یکی از چالشهای برجسته در طراحی کنترلکننده برای سیستمهای LPV تضمین پایداری سیستم در شرایطی است که پارامترهای دینامیکی آن به طور مداوم تغییر میکنند. به منظور دستیابی به این هدف، پایداری به وسیلۀ معیار لیاپانوف بررسی میشود و LMIها به عنوان ابزارهای اصلی و اثباتپذیر در تحلیل پایداری سیستمهای LPV به کار میروند ]17[. در این سیستمها، تابع لیاپانوف که میتواند به پارامتر وابسته باشد، به صورت تابعی از متغیرهای سیستم تعریف میشود که رابطۀ پایداری سیستم را در مواجهه با تغییرات دینامیک پارامترها مشخص میکند ]18[. در ]19[، یک رویکرد مدل مرجع ارائه شده است که در آن مدل خطای غیرخطی حاصل به یک مدل شبهخطی تغییرپذیر با پارامتر برای طراحی یک کنترلکنندۀ qLPV با استفاده از تکنیکهای مبتنی بر LMI آورده شده است.
3-2- مدلسازی رباتهای سیار غیرهولونومیک
رباتهای سیار چرخدار در بسیاری از زمینههای صنعتی و پژوهشی کاربردی گسترده دارند ]20، 21[. این رباتها به دلیل ویژگیهای خاص ساختاری خود، با قیدهای غیرهولونومیک مواجه هستند ]22[ که این قیدها حرکت آنها را محدود میکنند و در نتیجه، به طور طبیعی، آزادی حرکت آنها در محیطهای مختلف را کاهش میدهند ]23[. مدل سینماتیکی برای رباتهای غیرهولونومیک که نشاندهندۀ روابط پیچیده بین ورودیهای کنترلی و حرکت ربات است، به شکل زیر تعریف میشود:
(1) {█(x ̇= v cos(θ)@y ̇= v sin(θ)@θ ̇= ω)┤
که در آن (x,y) موقعیت ربات، θ زاویۀ حرکت، v سرعت خطی (v= 1/2(v_r+v_l))، v_r سرعت چرخ سمت راست، v_l سرعت چرخ سمت چپ وω سرعت زاویهای ربات (ω=1/2l(v_r-v_l)) است. این مدل به عنوان پایهای برای طراحی کنترلکنندۀ LPV در این مقاله استفاده شده است. شکل (1) شماتیک یک ربات سیار را نشان میدهد.
شکل (1): شماتیک ربات سیار
4-2- شناسایی و جبران خطای محرکها با استفاده از فیلتر کالمن توسعهیافته
در بسیاری از سیستمهای رباتیکی، خطاهای ناشی از مشکلات محرکها، همچون خرابیهای مکانیکی، سایش اجزا یا تغییرات در شرایط محیطی، میتوانند به طرزی جالب توجه بر کارایی و دقت سیستم تأثیر منفی داشته باشند ]24[. استفاده از فیلتر کالمن توسعهیافته به عنوان یک رویکرد کارآمد و رایج برای شناسایی و تخمین دقیق این خطاها شناخته میشود. این روش با اعمال تکنیکهای پیشرفته برای تخمین وضعیت سیستم در حضور خطا، امکان جبران خطاها و بهبود عملکرد را فراهم میآورد ]25، 26[. مدل دینامیکی ربات متحرک در شرایطی که خطای محرک در آن دخیل است، به صورت زیر تعریف میشود:
(2) {█(ω_L= (1 - δ_L ) ω_L^0@ω_R= (1 - δ_R ) ω_R^0 )┤
که در آن (δ_L,δ_R) ضرایب خطای محرک چپ و راست هستند و ω_L (t) و ω_R (t) بهترتیب سرعتهای زاویهای چرخهای چپ و راست را نشان میدهند. برای تخمین این مقادیر، فیلتر کالمن توسعهیافته استفاده میشود که با استفاده از یک مدل غیرخطی زمانگسسته، تخمینی دقیق از پارامترهای خطا ارائه میدهد. این روش در این مقاله برای بهبود عملکرد کنترلکنندۀ LPV پیشنهاد شده است.
در ]27[، دو روش برای تخمین موقعیت ربات پیشنهاد شدهاند و سپس یک فیلتر کالمن توسعهیافته برای ردیابی خروجیها و تشخیص هرگونه انحراف بزرگ غیرمنتظره استفاده میشود. در ]28[، خطاهای سیستم در رباتهای متحرک بر اساس مقایسۀ بین شبیهسازی و رفتار واقعی ربات تشخیص داده میشوند. یک الگوریتم تشخیص خطا برای سنسورها و محرکهای یک ربات متحرک محرک دیفرانسیلی در ]29[ ارائه شده است که بر اساس ماژولهای موقعیتیابی چندگانه است. یک رویکرد مبتنی بر مدل در ]30[ برای تشخیص خطاهای محرک با استفاده از رویکرد مبتنی بر ناظر استفاده میشود. هزینۀ زیاد به دلیل استفاده از سنسورهای موقعیتیابی اضافی عیب اصلی رویکرد افزونگی سختافزاری پیشنهادی است.
در این قسمت، فرایند تشخیص و شناسایی خطاها شرح داده میشود [31]. در این کار، از دست دادن خطای اثربخشی (LOE ) در محرکهای ربات متحرک در نظر گرفته میشود که میتواند به صورت زیر مدل شود:
ω_Rf (t)=(1-K_R)ω_R (t)
ω_Lf (t)=(1-K_L)ω_L (t) (3)
که در آن ω_Lf (t) و ω_Rf (t) بهترتیب سرعتهای زاویهای معیوب چرخهای چپ و راست را نشان میدهند و K_R و K_L بهترتیب از دست دادن مقادیر خطای مؤثر برای چرخهای راست و چپ را نشان میدهند. برای مثال، از دست دادن ۳۰ درصد از میزان اثر، محرک سمت چپ را میتوان با ضریب K_L=0.3 مدل کرد.
در ادامه، برای استفاده از این دو پارامتر باید بتوان آنها را تخمین زد که در این مقاله از فیلتر کالمن توسعهیافته استفاده شده است. برای دستیابی به این هدف، مدل غیرخطی زمانگسستۀ تکمیلی ربات متحرک با گسستهسازی رابطۀ (1) به روش اویلر و با در نظر گرفتن زمان نمونهبرداری به اندارۀ T، به صورت زیر ارائه شده است:
(4) [■(x(k+1)@y(k+1)@■(θ(k+1)@K_L (k+1)@K_R (k+1)))]=[■(x(k)+Tv(k)cos〖θ(k)〗@y(k)+Tv(k)sin〖θ(k)〗@■(θ(k)+Tω(k)@K_L (k)@K_R (k)))]
بر اساس معادلات (3)، میتوان (4) را به صورت زیر بازنویسی کرد:
(5) X(k+1)=f(X(k),u(k))
=[■(f_1 (X(k),u(k))@f_2 (X(k),u(k))@■(f_3 (X(k),u(k))@K_L (k)@K_R (k)))]
که در آن:
X(k)=[■(x(k)&y(k)&■(θ(k)&■(K_L (k)&K_R (k))))]^T
u(k)=[■(ω_L (k)&ω_R (k))]^T
f_1 (X(k),u(k))=x(k)+Tr/2((1-K_R (k)) ω_R (k)+(1-K_L (k)) ω_L (k))cos〖θ(k)〗
f_2 (X(k),u(k))=y(k)+Tr/2((1-K_R (k)) ω_R (k)+(1-K_L (k)) ω_L (k))sin〖θ(k)〗
f_3 (X(k),u(k))=θ(k)+Tr/d (1-K_R (k)) ω_R (k)+(1-K_L (k)) ω_L (k)
و T زمان نمونهبرداری است.
در این قسمت، از یک فیلتر کالمن توسعهیافته برای تخمین حالتهای X(k) به کمک سیستم غیرخطی زمان گسسته در رابطۀ (4) استفاده شده است که منجر به تخمین مقادیر خطای K_R و K_L میشود. برای رسیدن به این هدف، مدل غیرخطی زمان-گسستۀ زیر در نظر گرفته میشود، یعنی:
(6) X(k)=f(X(k-1),u(k-1))+w(k-1)
z(k)=h(X(k))+v(k)
که در آن X(k)∈R^(n_x ) بردار حالت، u∈R^(n_u ) بردار ورودی، z(k)∈R^(n_z ) بردار خروجی و h یک تابع غیرخطی با ابعاد مناسب است. w(k-1) و v(k) نویزهای پروسه و مشاهده هستند که بهترتیب، مستقل از فرایند تصادفی گاوسی سفید با متوسط صفر با ماتریسهای کوواریانس R ̅ و(Q ) ̅ فرض میشوند.
EKF توسط مرحلۀ پیشبینی با توجه به موارد زیر مشخص میشود:
(7) X ̂_(k|k-1)=f(X_(k-1),u_(k-1) )
P_(k|k-1)=F_(k-1) P_(k-1) F_(k-1)^T+Q ̅_(k-1)
در معادلات بالا، X ̂_(k-1) بردار حالت تخمینزدهشده در نمونۀ k-1 و P_(k-1) ماتریس کوواریانس در نمونۀ k-1 است و X ̂_(k|k-1) متغیر حالت در نمونۀ k را بر اساس اطلاعات سیستم تا نمونۀ k-1 پیشبینی میکند و Pk|k-1 پیشبینی ماتریس کوواریانس در نمونۀ kبر اساس اطلاعات سیستم تا نمونۀ k-1 است. ماتریس F_k دینامیک خطیشده است که به وسیلۀ خطیسازی تیلور به صورت زیر به دست میآید:
(8) F_k=∂f/∂X |_(X_(k-1) )
H_(k ) همچنین دینامیک اندازهگیری خطی استفادهشده در فاز بهروزرسانی است. این ماتریس نیز به صورت زیر محاسبه میشود:
(9) H_k=∂h/∂X |_(X ̂_(k|k-1) )
فاز بهروزیافتۀ فیلتر کالمن توسعهیافته با استفاده از رابطههای زیر انجام میشود:
(10) K_k=P_(k-1) H_k^T (H_k P_(k|k-1) H_k^T+R ̅_k)
X ̂_k=X ̂_(k|k-1)+K_k [z_k-h(X ̂_(k|k-1))]
P_k=(1-K_k H_k)P_(k|k-1)
5-2- کنترل ربات سیار دارای اشباع ورودی
یکی از مشکلات اساسی و رایج در طراحی سیستمهای کنترلی برای رباتهای متحرک محدودیتهای ناشی از اشباع ورودی هستند که به ویژه در محرکهای الکتریکی و مکانیکی مشاهده میشوند. زمانی که فرمانهای کنترلی از حدود مجاز خود فراتر میروند، این امر میتواند موجب بروز ناپایداری در سیستم یا کاهش دقت ردیابی مسیر شود ]32[. به همین دلیل، مدلسازی دقیق اشباع ورودی امری ضروری است تا از بروز چنین مشکلاتی جلوگیری شود ]33[. مسئلۀ در نظر گرفتن اشباع در حین طراحی کنترلکننده در سیستمهای کنترلی یکی از قیود مهم و با ارزش عملی زیاد است که باعث پیچیدگی در طراحی نیز میشود. با لحاظکردن این نکته در حین طراحی، از ناپایداری سیستم عملی به دلیل عدم اعمال سیگنال کنترلی مدنظر کنترلکننده به خاطر وجود اشباع جلوگیری میشود. در ادامه، مدل ریاضی مرتبط با اشباع ورودی به صورت زیر تعریف میشود:
(11) u_sat (u) = {█(u_max , u>u_max@u ,u_min≤u≤u_max@u_min ,u
که در آن u_max و u_min بهترتیب سطوح اشباع مجاز بالا و پایین هستند. برای در نظر گرفتن این محدودیت در طراحی کنترلکنندۀ LPV، رویکردهای مبتنی بر نامساویهای ماتریسی خطی به کار گرفته شدهاند.
6-2- نواحی پایداری تعمیمیافته
در برخی از کاربردها در مهندسی کنترل ممکن است نواحی پایداری خاصی برای ما اهمیت بیشتری داشته باشند. برای این منظور و با هدف بهبود پاسخ گذرای سیستم، در زمان طراحی، محل قطبهای سیستم تعیین میشود.
سیستم خطی زیر:
ẋ = Ax (12)
به طور مجانبی پایدار است، اگر و فقط اگر همۀ مقادیر ویژۀ A در -ℂ قرار گیرند. ناحیۀ پایداری به عنوان یک زیرمجموعۀ C_stab ⊆ C به صورت زیر تعریف میشود ]34[:
λ∈C_stab ⇒ 〖λ ̅∈C〗_stab (13)
C_stab is convex. (14)
ناحیۀ دلخواه در این مقاله با استفاده از سه تعریف زیر به دست آورده شده است:
C_(stab 1)= {s∈ C | Re{s}<-α} (15)
C_(stab 2)= {s∈ C | |s|<r} (16)
C_(stab 3)= {s∈ C | Re{s} tanθ<-|Im(s)|} (17)
رابطۀ (15) ناحیه را در سمت چپ خط عمودی δ = -α قرار میدهد. رابطۀ (16) یک ناحیۀ دایرهای با شعاع r در مرکز مبدأ را تعریف میکند. رابطۀ (17) یک بخش مخروطی را با رأس آن در مبدأ صفحۀ مختلط و زاویۀ داخلی θ تعریف میکند.
قضیۀ (2-1) ]34[: تمامی مقادیر ویژۀ A∈R^(n×n) در ناحیۀ LMI زیر قرار میگیرند:
{s∈ C | (■(I@sI))^* (■(Q&S@S^T&R))(■(I@sI))}<0 (18)
اگر و فقط اگر P>0 وجود داشته باشد، به گونهای که:
{ (■(I@A⊗I))^* (■(P⊗Q&P⊗S@P⊗S^T&P⊗R))(■(I@A⊗I))}<0 (19)
که در آن عملگر (⊗) نشاندهندۀ ضرب کرونکر است که به صورت زیر تعریف میشود:
(20) A⊗B=[■(A_11 B&⋯&A_1n B@⋮&⋱&⋮@A_m1 B&…&A_nm B)]
طراحی ردیاب LPV در ربات سیار دارای خطای محرک و اشباع ورودی
1-3- استخراج مدل LPV ربات سیار
در این بخش، یک مدل مرجع برای به دست آوردن مدل سینماتیکی توسعه داده خواهد شد. این مدل به صورت تفاوت بین اندازهگیریهای مقادیر واقعی (x, y, θ) و مقادیر مطلوب (xr, yr, θr) تعریف میشود که خطای موقعیت نامیده میشوند. شکل (2) ربات مرجع یک ربات فرضی است که به طور ایدهآل مسیر مرجع را دنبال میکند. در مقابل، ربات واقعی (در مقایسه با ربات مرجع)، خطا هنگام ردیابی مسیر مرجع است. بنابراین، الگوریتم کنترل باید طوری طراحی شود که ربات، مسیر مرجع را به طور دقیق دنبال کند.
شکل (2): نمایی از خطای ردیابی ربات
بنابراین، با استفاده از معادلات سینماتیک ربات (1)، خطای وضعیت فعلی ربات متحرک نسبت به حالت مورد انتظار در بدنۀ ربات را میتوان به صورت زیر به دست آورد:
(21) [■(e_1@e_2@e_3 )]=[■(cosθ&sinθ&0@〖-sin〗θ&cosθ&0@0&0&1)][■(x_r-x@y_r-y@θ_r-θ)]
از معادلۀ بالا مشتق گرفته شده و با استفاده از مدل سینماتیکی ارائهشده در بخش دوم در رابطۀ (2)، دینامیک خطای ربات محاسبه میشود:
(22) [■(e ̇_1@e ̇_2@e ̇_3 )]=[■(cos〖e_3 〗&0@sin〖e_3 〗&0@0&1)][■(v_r@ω_r )]+[■(-1&e_2@0&-e_1@0&-1)][■(v@ω)]
معمولا از یک کنترل سرعت ورودی استفاده میشود که پایداری (1) را تضمین میکند [35]:
(23) {■(v=v_r cos〖e_3 〗-v_v@ω=ω_r-v_ω )┤
که در آن 〖 v〗_vو v_ωسیگنالهای کنترلی هستند که باید توسط کنترلکننده محاسبه شوند. و 〖 v〗_r 〖,cos〗〖e_3 〗و ω_rسیگنالهایی هستند که توسط کنترلکنندۀ پیشخور تولید میشوند. استفاده از ورودیهای (13) در (14) مدل غیرخطی زیر را نتیجه میدهد:
(24) [■(e ̇_1@e ̇_2@e ̇_3 )]=[■(0&ω_r (t)&0@-ω_r (t)&0&0@0&0&0)][■(e_1@e_2@e_3 )]+[■(0@sin〖(e_3)〗@0)] v_r+[■(1&0@0&0@0&1)][■(v_v@v_ω )]
در نهایت، معادلۀ (24) را میتوان به صورت زیر نیز نوشت:
(25)
e ̅ ̇=[■(0&ω_r (t)&0@-ω_r (t)&0&v_r (t)(sin(e_3 )/e_3 )@0&0&0)] e ̅+[■(1&0@0&0@0&1)].[■(v_v@v_ω )]
که در آن e ̅ ̇ بردار دینامیک خطاست، 〖 v〗_rو ω_rسرعتهای خطی و زاویهای مرجع هستند، و v_v وv_ω سرعتهای خطی و زاویهای ربات متحرک به عنوان ورودی سیستم هستند. همانطور که در معادلۀ (25) دیده میشود، معادلات خطیشده دارای پارامترهایی هستند که با زمان تغییر میکنند. این پارامترها با استفاده از روش LPV تعریف شدهاند. سیستم LPV به صورت زیر توصیف میشود:
(26) e ̅ ̇=A(δ) e ̅+Bu
=[■(0&δ_1&0@-δ_1&0&δ_2@0&0&0)] e ̅+[■(1&0@0&0@0&1)].[■(v_v@v_ω )]
که δ_2∈[-0.2*v_rmax ,v_rmax ] و δ_1∈[ω_rmin,ω_rmax ] است.
2-3- طراحی کنترلکنندۀ LPV برای ربات بدون در نظر گرفتن اشباع ورودی
فرض کنید سیستم LPV در رابطۀ (26) به صورت زیر در نظر گرفته شده است:
(27) e ̅ ̇=A(δ) e ̅+Bu
که در آن e ̅∈R^(n_e ) بردار خطا، u∈R^(n_u ) بردار ورودی است و δ∈∆ بردار پارامترهای متغیر است که در آن:
(28) ∆=Conv(∆_0 )
∆_0={δ∈R^(n_δ ) |δ_k∈{▁(δ_k ),¯(δ_k )}}
از مدل (27) برای طراحی کنترلکنندۀ فیدبک حالت برای مدل پیشنهادی استفاده شده است. با استفاده از فیدبک حالت زیر داریم:
(29) u=Ke ̅
هدف یافتن یک بهرۀ کنترلکنندۀ K است؛ به نحوی که پایداری و عملکرد را در ناحیۀ چندوجهی تضمین کند. این هدف را میتوان با استفاده از پایداری لیاپانوف، نابرابریهای ماتریسی خطی (LMI) و ابزارهای بهینهسازی محدب به دست آورد.
با توجه به سیستم حلقهبستۀ e ̅ ̇=(A(δ)-BK) e ̅، شرط کافی برای تضمین پایداری سیستم وجود یک تابع درجهدوم به فرم زیر است:
(30) V(e ̅ )=e ̅^T We ̅>0
به نحوی که رابطۀ زیر برای مشتق آن برقرار باشد:
(31) V ̇(e ̅ )=e ̅^T [(A(δ)-BK)^T W+W(A(δ)-BK)] e ̅<0
که در آن W یک ماتریس مثبت معین متقارن است [36].
نامعادلۀ (31) را میتوان به شکل LMI به صورت زیر بازنویسی کرد:
(32) W>0
A(δ)W+W〖A(δ)〗^T-BX-X^T B^T <0
پس از حل LMIهای بالا، بهرۀ کنترلی از رابطۀ K =XW^(-1) به دست میآید که منجر به کنترلکنندۀ مطلوب میشود و قطبهای حلقهبسته را در نیمصفحۀ سمت چپ قرار میدهد. با این حال، حل این نامعادله فقط منجر به قرارگیری قطبهای حلقهبسته در سمت چپ محور موهومی خواهد شد. محدودکردن ناحیهای که در آن قطبهای حلقهبسته با استفاده از بخش (2-6) قرار داده میشوند، ترجیح داده میشود؛ به طوری که برخی از معیارهای عملکرد حاصل شوند.
3-3- طراحی کنترلکنندۀ LPV برای ربات با در نظر گرفتن اشباع ورودی و ناحیۀ پایداری LMI
همانطور که در بخش (2-5) و به کمک ]37[ اشاره شد، رابطۀ زیر برقرار است:
(33) e ̅ ̇=A(δ) e ̅+B{sat(u(t))}
که در آن e ̅ و u(t) بهترتیب بردار خطا و ورودی کنترل هستند. علاوه بر این، sat(.) نشاندهندۀ عملگر اشباع است که به صورت زیر تعریف میشود:
(34) [sat(u)]_i≜{█([u]_i; [u]_i< μ @μ; [u]_i≥ μ@-μ; [u]_i≤ -μ)┤
که در آن [u]_i به عنوان iامین عنصر از u و μ نیز سطح اشباع تعریف میشود.
لم (3-1) ]38[: با در نظر گرفتن u,v∈R^mبه صورت زیر:
(35) u=[u_1,u_2,…,u_m ]^T
v=[v_1,v_2,…,v_m ]^T
و با فرض |v_i | ≤μ برای همۀ i∈[1,m]، رابطۀ زیر برقرار است:
(36) sat(u)∈co{D_i u+D ̅_i v}
که در آن co نشاندهندۀ ناحیۀ محدب و D مجموعهای از ماتریسهای قطری m×m است که عناصر قطری آن 0 یا 1 هستند. سیستم بالا بهازای m=2 برقرار است؛ به طوری که:
D={[■(0&0@0&0)],[■(0&0@0&1)],[■(1&0@0&0)],[■(1&0@0&1)]} ,
D ̅_i=I-D_i
طبق لم (3-1)، میتوان رابطۀ زیر را به دست آورد:
(37) sat(u)= ∑_(j=1)^(2^m)▒〖ξ_j {D_j u+D ̅_j v} 〗
∑_(j=1)^(2^m)▒〖ξ_j=1 〗 ,ξ_j≥0
تعریف (3-1) ]37[: برای یک اسکالر مثبت ρ و یک ماتریس متقارن P>0 و با تعریف V_x=X^T (t)PX(t) و Ω(P,ρ)=〖{X∈R^n,X〗^T PX<ρ}، مجموعۀ Ω(P,ρ) مجموعهای نامتغیر است اگر رابطۀ V ̇_x<0 برای همۀ X∈Ω(P,ρ) برقرار باشد.
از مدل (33) برای طراحی کنترلکنندۀ فیدبک حالت برای مدل پیشنهادی استفاده شده است. از رابطۀ زیر برای فیدبک حالت استفاده شده است:
(38) u=Ke ̅
قضیۀ (3-1): سیستمe ̅ ̇=A(δ) e ̅+Bsat(Ke ̅ ) به صورت مجانبی پایدار است و تمامی مقادیر ویژۀ A(δ)∈R^(n×n) در ناحیۀ LMI زیر قرار میگیرند، اگر:
(39) ∃ P_1>0 & H_1,K_1
s.t. (A(δ) P_1+B ̂_η )^T+(A(δ) P_1+B ̂_η )<0 ,
B ̂_η=B∑_(j=1)^4▒〖ξ_j {D_j K_1+D ̅_j H_1 } 〗
(40) [■(U_sat^T U_sat&H_1@H_1^T&P_1 )]>0
U_sat =[■(μ_1&μ_2 )]^T
|V|<μ_1 ,|ω|<μ_2
(41) A(δ) P_1+P_1 〖A(δ)〗^T-B ̂_η
-B ̂_η^T+2αP_1 <0
(42) [■(-rP_1&P_1 〖A(δ)〗^T-B ̂_η^T@A(δ) P_1-B ̂_η&-rP_1 )]<0
(43) [■(sinφ A ̅_1&cosφ A ̅_2@cosφ A ̅_3&sinφ A ̅_4 )]<0
A ̅_1=A(δ) P_1+P_1 〖A(δ)〗^T-B ̂_η-B ̂_η^T
A ̅_2=A(δ) P_1-P_1 〖A(δ)〗^T-B ̂_η+B ̂_η^T
A ̅_3=A(δ) P_1-P_1 〖A(δ)〗^T+B ̂_η-B ̂_η^T
A ̅_4=A(δ) P_1+P_1 〖A(δ)〗^T-B ̂_η-B ̂_η^T
اثبات: هدف، یافتن بهرۀ کنترلکنندۀ K است که پایداری و عملکرد را در ناحیۀ چندوجهی تضمین کند. برای این امر، میتوان با استفاده از روابط زیر به این مهم دست یافت:
ابتدا تابع لیاپانوف زیر در نظر گرفته شده است:
(44) V(e ̅ )=e ̅^T Pe ̅>0
مشتق زمانی این تابع به صورت زیر به دست میآید:
(45) V ̇(e ̅ )=e ̅ ̇^T Pe ̅+e ̅^T Pe ̅ ̇<0
(46) V ̇(e ̅ )=[A(δ) e ̅+Bsat(Ke ̅ )]^T Pe ̅+e ̅^T P[A(δ) e ̅+Bsat(Ke ̅ )]
با استفاده از رابطۀ (37) و v(t)=He ̅ ,|h_i e ̅ |≤μ ، تابع اشباع به صورت زیر بازنویسی میشود:
(47) 𝑠𝑎𝑡(𝐾e ̅)= ∑_(j=1)^4▒〖ξ_j {D_j K+(D_j ) ̅ 〗 H}e ̅
∑_(j=1)^4▒〖ξ_j=1 〗 ,ξ_j≥0
بنابراین، مشتق تابع لیاپانوف به صورت زیر بازنویسی میشود:
(48) V ̇(e ̅ )=〖e ̅^T [A(δ)+B∑_(j=1)^4▒〖ξ_j {D_j K+ ( D) ̅_j H} 〗]〗^T Pe ̅+e ̅^T P[A(δ)+ B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗] e ̅<0
بنابراین، شرط کافی برای منفی معین شدن این تابع رابطۀ زیر است:
(49) ∑_(j=1)^4▒〖ξ_j ([A(δ)+B{D_j K+ D ̅_j H}]^T P+P[A(δ)+B{D_j K+ D ̅_j z}]) 〗<0
رابطۀ (49) خطی نیست و با ضرب عبارت P_1=P^(-1) در طرفین نامعادلۀ بالا و با تعریف K_1=KP_1و H_1=HP_1 و B ̂_η=B∑_(j=1)^4▒〖ξ_j {D_j K_1+D ̅_j H_1 } 〗 آنگاه رابطۀ (39) اثبات میشود.
با توجه به تعریف (3-1)، Ω(P_1,ρ) یک مجموعۀ نامتغیر از (44) برای همۀ e ̅∈Ω(P,ρ) است. اشباع در صورتی رخ نمیدهد که مجموعۀ نامتغیر درون L|H| باشد، یعنی Ω(P_1,ρ)⊂L|H|. با فرض ρ=1، رابطۀ زیر برقرار است:
(50) H^T PH<U_sat^T U_sat
〖⇒U〗_sat^T U_sat-H^T PH>0
با استفاده از لم شور میتوان به رابطۀ زیر رسید:
(51) [■(U_sat^T U_sat&H@H^T&P)]>0
با ضرب ماتریس[■(I&0@0&P^(-1) )] از سمت چپ و راست ماتریس بالا، رابطۀ (32) اثبات میشود.
رابطۀ (41) ناحیه پایداری را در سمت چپ خط عمودی δ = -α قرار میدهد که برای اثبات آن از قضیۀ (2-1) استفاده میشود:
(52) Q=2α ,S=1 ,R=0
با جایگذاری مقادیر بالا، روابط زیر حاصل میشوند:
(53) (■(I@A(δ)-B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗))^* (■(2αP&P@P&0))
(■(I@A(δ)-B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗))<0
(54) 2αP+(A(δ)-B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗)^T P+P(A(δ)-B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗)<0
با ضرب رابطۀ (54) از سمت چپ و راست در P_1=P^(-1)، رابطۀ (41) اثبات میشود.
رابطۀ (42) یک ناحیۀ دایرهای با شعاع r در مرکز مبدأ را تعریف میکند که برای اثبات آن از قضیۀ (2-1) استفاده میشود:
(55) Q=-r^2 ,S=0 ,R=1
B_1= B∑_(j=1)^4▒〖ξ_j {D_j K+D ̅_j H} 〗
(56) { (■(I@A(δ)-B_1 ))^* (■(-r^2 P&0@0&P))(■(I@A(δ)-B_1 ))}<0
(57) -r^2 P+(A(δ)-B_1 )^T P(A(δ)-B_1 )<0
در رابطۀ (57) از سمت چپ و راست P_1=P^(-1) ضرب و با استفاده از لم شور رابطۀ (33) اثبات میشود. رابطۀ (43) یک بخش مخروطی را با رأس آن در مبدأ صفحۀ مختلط و زاویۀ داخلی 2φ تعریف میکند:
(58) Q=[■(0&0@0&0)] ,S=[■(sinφ&cosφ@-cosφ&sinφ )] , R=[■(0&0@0&0)]
با جایگذاری مقادیر بالا در قضیۀ (2-1)، رابطۀ (43) نیز اثبات خواهد شد.
4-3- کنترل خطای محرک
همانطور که در بخش (2-5) شرح داده شد، مقادیر برآوردشدۀ K ̂_L و K ̂_R از الگوریتم EKF در این بخش برای حذف اثرات خطا بر روی عملکرد ردیابی مسیر رباتها استفاده خواهند شد.
برای اینکه کنترلکنندۀ طراحیشده مقاوم در برابر خطا باشد، ابتدا کنترلکنندۀ طراحیشده در بخش (2-4) در نظر گرفته میشود و سپس خروجی کنترلکننده به صورت زیر اصلاح میشود:
(59) ω_R (t)=1/(1-K ̂_R ) (v(t)/r+dω(t)/2)
(60) ω_L (t)=1/(1-K ̂_L ) (v(t)/r-dω(t)/2)
تأثیر خطای LOE با طرح بالا لغو میشود، با فرض اینکه ω_R (t)و ω_L (t) از حداکثر سرعت مجاز چرخ ربات تجاوز نکنند. بلوک دیاگرام در شکل (3) نشان داده شده است.
نکتۀ 1: در این بخش، ابتدا دینامیک سیستم خطا با رابطۀ (18) به فرم LPV درآورده شده است. سپس در رابطۀ (32) یک LMI بیان شده است که با استفاده از آن میتوان کنترلکننده را برای سیستم بدون در نظر گرفتن اشباع ورودی طراحی کرد. در ادامه، طبق رابطۀ (33) و با ارائۀ قضیۀ (3-1)، رابطۀ (43) به دست آمد که میتوان با استفاده از آن یک کنترلکننده برای حالت با در نظر گرفتن اشباع ورودی طراحی کرد. در نهایت، برای مقابله با اثر خطای محرک، روابط کنترلی (59) و (60) ارائه شدهاند.
شکل (3): بلوک دیاگرام کنترل خطای مبتنی بر LPV رباتهای متحرک
نکتۀ 2: در مقایسه با روش ارائهشده در [37]، روش ارائهشده در این مقاله نسبت به خطای محرک مقاوم است. همچنین، به عنوان مزیت دیگر، در طراحی کنترلکننده در این مقاله، نواحی پایداری LMI در نظر گرفته شده است که باعث بهبود پاسخ حالت گذرا میشود. همچنین، در مقایسه با [31]، از روش مدلسازی LPV استفاده شده است که باعث میشود مقاومت روش در برابر تغییرات پارامترها بیشتر باشد. در بخش بعد، با استفاده از شبیهسازی، مقایسهای دقیقتر با این روش انجام خواهد شد.
نتایج شبیهسازی
در این بخش، نتایج شبیهسازی ردیابی مسیر ربات سیار غیرهولونومیک با استفاده از کنترلکنندههای LPV ارائه میشود که در بخش قبل طراحی شدهاند.
1-4- طراحی کنترلکنندۀ LPV برای ربات با در نظر گرفتن اشباع ورودی
مدل LPV استخراجشده برای سیستم غیرخطی ربات متحرک، که در بخش سوم به آن اشاره شد، به صورت زیر تعریف میشود:
(61) e ̅ ̇=A(δ) e ̅+Bu
=[■(0&δ_1&0@-δ_1&0&δ_2@0&0&0)] e ̅+[■(1&0@0&0@0&1)].[■(v_v@v_ω )]
بازۀ تغییرات درایههای ماتریس بالا به صورت زیر است:
(62) -1.4≤δ_1≤1.4 ,-0.2≤δ_2≤0.5
ناحیۀ پایداری دلخواه که برای بهبود کیفیت پاسخ حالت گذرا در بخش سوم معرفی شد، مقادیر شعاع دایره و زاویۀ داخلی مخروطی و همچنین ناحیۀ سمت چپ خط عمودی α به صورت زیر در نظر گرفته شدهاند:
(63) α=-0.7 ,θ=80 ,r=3
در ادامه، با در نظر گرفتن روابط (61) تا (63) و بخش (3-2)، کنترلکنندۀ زیر طراحی میشود:
(64) K =[■(4.8&0.8&0.4@1.8&9.3&7.4)]+δ_1 [■(-0.5&-5.4&-1.2@-7.7&4.6&0.64)]+δ_2 [■(1. 3.5 0.9@-2.7 50 13.8)]
برای تخمین خطای محرک هم ماتریسهای کوواریانس نویز پروسه به صورت زیر در نظر گرفته شدهاند:
(65) Q ̅=6×10^(-6) I_6 ,R ̅=0.001I_3
در ادامه، با استفاده از محیط شبیهسازی سیمولینک، مدل اصلی ربات سیار غیرهولونومیک اجرا و سپس با کنترلکنندۀ (64) مسیر ردیابی شده است.
در این راستا، یک بار فرض شده است محرک سمت راست 30 درصد و یک بار هم سمت چپ 80 درصد کارایی خود را از دست دهد؛ بنابراین، نتایج را میتوان در شکلهای (4) و (5) مشاهده کرد:
شکل (4): ردیابی مسیر مرجع، زمانی که در t=10s محرک سمت راست 30 درصد کارایی خود را از دست دهد.
شکل (5): ردیابی مرجع، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست دهد.
در هر دو حالت، کنترلکنندۀ پیشنهادی بهخوبی توانسته است ربات را به مسیر مرجع بازگرداند، اگرچه در خطای بالاتر، دنبالکردن مسیر مرجع با چالش بیشتری مواجه است.
در شکل (6)، سیگنال کنترلی بهاشباعرفته نیز نشان داده شده است:
شکل (6): سیگنال کنترلی اشباع شده، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست دهد.
در ادامه، اثر ناحیۀ پایداری دلخواه که در طراحی کنترلکننده در قضیۀ (3-1) در نظر گرفته شده بود، با حالت بدون ناحیۀ پایداری دلخواه مقایسه شده است. در این راستا، پاسخ در حالتی که محرک سمت چپ 30 درصد کارایی خود را از دست داده است، با هم مقایسه میشود:
شکل (7): نمودار خطای ردیابی زمانی بدون در نظر گرفتن ناحیۀ پایداری
شکل (8): نمودار خطای ردیابی زمانی با در نظر گرفتن ناحیۀ پایداری
با توجه به شکلهای (7) و (8)، واضح است میزان فراجهشها و فروجهشها و همچنین زمان نشست کاهش پیدا کرده است.
در ادامه، تأثیر در نظر گرفتن/نگرفتن شرط اشباع بر روی سیگنال کنترلی در حین طراحی ارائه میشود. در شکل (9)، ردیابی مسیر مرجع با استفاده از کنترلر LPV بدون شرط اشباع، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست بدهد، نشان داده شده است. در این شبیهسازی، فرض شده است در هنگام طراحی، شرط اشباع بر سیگنال کنترل در نظر گرفته نشده است، اما در حین شبیهسازی، محرکها به اشباع میروند. همانطور که از این شکل پیداست، ربات از مسیر خود خارج شده و به مسیر اصلی بازنگشته است. بنابراین، لزوم در نظر گرفتن شرط اشباع در حین طراحی در اینجا مشخص میشود.
شکل (9): ردیابی بدون در نظر گرفتن شرط اشباع، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست بدهد.
3-4- مقایسه با روش [31]
به دلیل اینکه در این مقاله از روش ارائهشده در ]31[ برای کنترل خطا استفاده شده است، در ادامه، نتایج مقالۀ حاضر با این مقاله مقایسه خواهد شد. در ادامه، از آنجا که مرجع ]31[ با مسیر دایرهای نتایج خود را بررسی کرده است، در اینجا نیز از مسیر دایرهای برای مقایسه استفاده خواهد شد و فرض شده است محرک سمت چپ 80 درصد کارایی خود را از دست دهد. نتایج را میتوان در شکلهای (10) و (11) مشاهده کرد.
با بررسی این نتایج، میتوان دریافت ردیابی مسیری که کنترلر LPV با در نظر گرفتن اشباع ورودی از خود به نمایش گذاشته است دارای عملکردی بهتر نسبت به ردیابی کنترلکنندۀ ]31[ است. گفتنی است، روش ارائهشده در مرجع ]31[ در ردیابی مسیر پروانهای ناموفق بوده است.
شکل (10): ردیابی مسیر مرجع با استفاده از کنترلکنندۀ ]31[، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست دهد.
شکل (11): ردیابی مسیر مرجع با کنترلکنندۀ LPV با در نظر گرفتن اشباع ورودی پیشنهادی، زمانی که در t=10s محرک سمت چپ 80 درصد کارایی خود را از دست دهد.
نتیجهگیری
در این مقاله، کنترلکنندۀ مبتنی بر مدل LPV برای رباتهای متحرک غیرهولونومیک طراحی شد. این کنترلکننده قابلیت جبران اثرات اشباع ورودی و خطاهای ناشی از آن را دارد. نتایج شبیهسازیها نشان داد این کنترلکننده به طرزی جالب توجه دقت ردیابی مسیر رباتها را بهبود بخشیده و توانسته است عملکرد سیستم را حتی در مواجهه با اشباع ورودی و اختلالات دیگر حفظ کند. با این حال، چالشهایی مانند پیچیدگیهای محاسباتی و نیاز به تنظیم دقیق پارامترها در محیطهای واقعی همچنان وجود دارند که به پژوهشهای بیشتر برای بهبود عملکرد این سیستمها در شرایط عملیاتی پیچیدهتر نیاز دارند.
| ||
| مراجع | ||
| ||
|
آمار تعداد مشاهده مقاله: 129 تعداد دریافت فایل اصل مقاله: 35 |
||