\documentclass[11pt,a4paper]{article}

%\usepackage{cite}
\usepackage[numbers]{natbib} 
\usepackage{float} 
%\usepackage[caption = false]{subfig} 
\usepackage{subfigure} 
\usepackage{caption}
\usepackage{multicol}
\usepackage{epsf,graphicx}
\usepackage{pdfsync}
\usepackage{graphicx}
\usepackage{pstricks}
\usepackage{color}
\usepackage{setspace}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{amsthm}
\usepackage[fleqn]{amsmath}
\usepackage{mathrsfs}
\usepackage{mathtools}

%\usepackage{multicol, latexsym, amsmath, amssymb} 
\usepackage{bm}

\usepackage{float} 
%**************
%for table
\usepackage{tabularx,booktabs,caption}
%*************

 
\usepackage[colorlinks,citecolor=blue]{hyperref}
\usepackage[top=25mm, bottom=25mm, left=20mm, right=20mm]{geometry}

\setlength{\columnwidth}{82mm}
\setlength{\columnsep}{6mm}

%\usepackage{bidipoem}

%\usepackage[Kashida]{xepersian}
\usepackage{xepersian}
\settextfont[Scale=1]{XB Zar}%{B Nazanin}%
\setlatintextfont[Scale=.9]{Times New Roman}
%\setdigitfont{XB Zar}
\defpersianfont\iranic[Scale=.8]{XB Zar Italic}
%\defpersianfont\Keyhan[Scale=1]{XB Kayhan Sayeh}

\defpersianfont\TitleBold[Scale=1]{XB Niloofar Bold}
\defpersianfont\AbstractBold[Scale=.92]{XB Niloofar Bold}
\defpersianfont\nastaliq[Scale=1.2]{IranNastaliq}
%\usepackage{persianpoem}


\newcommand{\textblue}[1]{{\addfontfeature{Color=0000FF}#1}}
\numberwithin{table}{section}
\onehalfspacing

\newcommand\femph[1]{\lr{''}#1\lr{``}}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclareMathOperator*\argmin{arg\,min}
\DeclareMathOperator*\argmax{arg\,max}
%\renewcommand{\algorithmicrequire}{\textbf{ورودی:}}
%\renewcommand{\algorithmicensure}{\textbf{خروجی:}}
\DefaultMathsDigits
\begin{document}
\SepMark{-}	

\title{\TitleBold{کنترل غیرخطی مقید ربات اسکارا با استفاده از کنترل‌کننده تطبیقی IQ-PD}}

\author{نجمه اسکندری$^{1}$، ...$^2$، ...$^3$ و ...$^4$\\
$^\dag$دانشگاه صنعتی اصفهان، دانشکده مهندسی برق و کامپیوتر\\
}

%\twocolumn[
\thispagestyle{empty}

%%%%%%%%% ABSTRACT%%%%%%%%%
%\begin{@twocolumnfalse}
\date{}
\pagestyle{empty}
\maketitle\thispagestyle{empty}
\begin{abstract}
\AbstractBold
{یکی از مهم‌ترین چالش‌های استفاده از ربات‌ها، حفظ ایمنی هنگام تعامل ربات با محیط اطراف خود به خصوص انسان‌ها است. در همین راستا، با ترکیب و طراحی دو کنترل‌کننده امپدانس و نامتغیر برای ربات اسکارا سعی شده است که تا حد امکان از آسیب به اپراتورها هنگام همکاری نزدیک با ربات و همچنین از آسیب به خود ربات و ایجاد هزینه‌های گزاف جلوگیری شود. به دلیل عدم قطعیت در مدل‌سازی دینامیکی ربات و شرایط کاری متغیر، روش‌های کلاسیک تنظیم ضرایب کنترل‌کننده‌های تناسبی-مشتق‌گیر، که به منظور ردیابی مسیر در بلوک دیاگرام پیشنهادی قرار دارند، ممکن است به خوبی جوابگو نباشد. بنابراین در این مقاله در عین حفظ ایمنی، برای برآورده شدن هدف نهایی کنترل که همان ردیابی مسیر مطلوب می‌باشد، ضرایب بهینه با استفاده از الگوریتم یادگیری کیو افزایشی به صورت برخط و حین انجام عملیات توسط ربات، محاسبه می‌شوند.  }



 \end{abstract}
 

 \hspace{9mm}\textbf{کلمات کلیدی}: بازوهای مکانیکی، کنترل امپدانس، کنترل نامتغیر، یادگیری کیو افزایشی.\\

%\end{@twocolumnfalse}]
\pagestyle{myheadings}
\begin{multicols}{2}\raggedcolumns
\section{مقدمه}

پیشرفت‌های اخیر در حوزه رباتیک باعث استفاده بیشتر از ربات‌ها در حوزه‌های کاربردی متفاوت شده است. افزایش توانایی‌های محاسباتی، پیشرفت ساختاری گیره‌های ماهر\LTRfootnote{Dextrous Gripper} و بازوهای مکانیکی\LTRfootnote{Manipulators} قابلیت انجام عملیات‌های پیچیده را در محیط‌های پویا و در ارتباط با انسان ممکن می‌سازد. علی‌رغم کاربرد گسترده بازوهای مکانیکی در صنعت، استفاده از این سیستم‌ها در محیط‌های انسانی به شدت محدود می‌باشد. با این وجود تلاش‌های بسیار زیادی برای پیوند انسان و ربات انجام شده است. از جمله مزایای همکاری، استفاده از توانایی ربات‌های نیرومند برای انجام وظایف سنگین با هدایت انسان است \cite{1}. بنابراین در تعامل انسان و ربات باید  عوامل دقیق تهدیدکننده‌ی ایمنی شناسایی و محدودیت‌های لازم برای جلوگیری از خطر اعمال شود \cite{2}. 

به منظور برقراری ایمنی باید تمهیدی اندیشیده شود که از برخورد ناخواسته ربات به محیط اطراف جلوگیری کند و یا در صورت برخورد شدت آسیب وارده به محیط و ربات به حداقل برسد. راهکارهای جلوگیری از آسیب به دو دسته پیش از برخورد\LTRfootnote{Pre-Collision} و حین برخورد\LTRfootnote{Post-Collision} تقسیم می‌شوند. راهکارهای پیش از برخورد از برخورد ناخواسته ربات به موانع محیطی جلوگیری می‌کنند \cite{3,4}. هدف در راهکارهای حین برخورد، کنترل نیرو و کاهش شدت آسیب در صورت برخورد است. کاهش شدت آسیب از دو جنبه تعامل امن از طریق برنامه‌ریزی و کنترل\LTRfootnote{Interaction Safety Through Planning and Control} \cite{7,6,5} و تعامل امن از طریق طراحی ساختار ربات\LTRfootnote{Interaction Safety Through Design} مورد بررسی قرار می‌گیرد \cite{8}.  

در گام اول برقراری ایمنی، برای جلوگیری از برخورد ناخواسته ربات به محیط اطراف و اعمال محدودیت روی مکان مجری نهایی یا به طور کلی اعمال قید روی سیگنال کنترل، خروجی و متغیرهای حالت، روش‌های کنترلی متعددی چون کنترل پیش‌بین مدل\LTRfootnote{Model Predective Control} در مرجع \cite{9}، روش گاورنر 
مرجع\LTRfootnote{Reference Governor} در \cite{10}، کنترل نامتغیر\LTRfootnote{Invariance Control} و روش میدان پتانسیل\LTRfootnote{Potential Field} ارائه شده‌اند \cite{11}. در روش کنترل پیش‌بین مدل که بیشتر برای اعمال قید روی سیگنال کنترلی به کار برده می‌شود، ابتدا رفتار آینده سیستم روی افق زمانی محدودی پیش‌بینی و سپس در هر لحظه، سیگنال ورودی آینده با کمینه‌سازی یک تابع هدف تحت قیود ورودی محاسبه می‌شود. در روش گاورنر مرجع با فیدبک از حالت‌های سیستم، پیش‌بینی و بهینه‌سازی تابع هدف موقعیت مرجع طوری تغییر داده می‌شود که هیچ محدودیتی نقض نشود. به دلیل بهینه‌سازی برخط و بار محاسباتی، ممکن است این روش در کاربردهای عملی به خوبی جوابگو نباشد و این در حالی است که در کنترل نامتغیر وجود جواب تحلیلی بهینه‌سازی محدب از بار محاسباتی می‌کاهد. روش میدان پتانسیل با در نظر گرفتن ربات در میدان پتانسیل مجازی فرض می‌کند که موقعیت مطلوب نیروی جاذبه و موانع محیطی نیروی دافعه دارند و برآیند این دو نیرو بر ربات عمل می‌کند.از معایب این روش احتمال گیر افتادن ربات در مینیمم‌های محلی می‌باشد.  

کنترل نامتغیر اولین بار در \cite{12} برای کنترل و پایدارسازی سیستم‌های غیرخطی غیر مینیمم فاز و تک ورودی-تک خروجی مطرح شد. وولف\LTRfootnote{Wolff} در سال 2004 کنترل‌ نامتغیر را به عنوان یک کنترل با قابلیت اضافه شدن\LTRfootnote{Add-On Control} به کنترل نامی که بر اساس هدف اصلی سیستم طراحی می‌شود، معرفی کرد \cite{13}. این کنترل‌کننده با تعریف تابع محدودیت که معیاری از فاصله ربات تا محدودیت است، سعی دارد از مجموعه‌ حالتهای سیستم که مقدار تابع محدودیت به ازای آنها همواره منفی است و مجموعه نامتغیر نامیده می‌شود، خارج نشود \cite{14}.   


 اگر رباتی که به طور مثال مشغول پاک کردن شیشه است صرفاً کنترل موقعیت شود، در صورتی که موقعیت مرجع به اشتباه پشت شیشه تنظیم گردد، ربات بدون توجه به محیط و آسیب‌های احتمالی مقدار نیرو را تا جایی افزایش می‌دهد که در موقعیت مرجع قرار گیرد و بدین ترتیب شیشه می‌شکند. بنابراین نیازمند روشی هستیم که هنگام تعامل با محیط و یا به عبارت دیگر تبادل نیرو، از آسیب‌های احتمالی جلوگیری شود. کنترل امپدانس که در این مقاله به عنوان کنترل‌کننده‌ی نیرو مورد استفاده قرار گرفته است و در دسته راهکارهای حین برخورد قرار می‌گیرد با مدل کردن رفتار دینامیکی ربات به عنوان یک سیستم جرم-فنر-دمپر ضمن کنترل موقعیت، نیروی وارده از طرف ربات به محیط را کنترل می‌کند. بدین ترتیب در صورتی که ربات برخوردی ناگهانی با محیط اطراف داشته باشد، در مقابل تغییر مسیر مرجع مقاومتی از خود نشان نمی‌دهد و با حداقل انحراف لازم از مسیر، شدت برخورد کاهش می‌‌یابد \cite{15}.

کنترل سختی\LTRfootnote{Stiffness-Control} که در مرجع \cite{16} بیان شده‌ است، از جمله روش‌های مطرح در حوزه کنترل نیرو می‌باشد. لزوم ثابت بودن موقعیت یا نیروی مطلوب در کنترل سختی یکی از معایب اصلی این کنترل‌کننده است. بنابراین در مسائل کنترلی که هدف در آنها ردیابی یک مسیر مطلوب متغیر با زمان است نمی‌تواند مورد استفاده قرار بگیرد و این در حالی است که از کنترل امپدانس، هم در مسائل کنترلی تنظیم\LTRfootnote{Regulation} و هم در مسائل ردیابی\LTRfootnote{Tracking} می‌توان استفاده نمود.




عملکرد کنترل‌کننده‌های PD که برای ردیابی مسیر مرجع اصلاح‌شده استفاده می‌شوند، به شدت به تنظیم مناسب ضرایب وابسته است. به دلایلی چون شرایط کاری متغیر، عدم قطعیت در مدل‌سازی دینامیکی و تغییر پارامترهای دینامیکی ربات (به علل متعددی مانند بلند کردن اجسام با جرم‌های متفاوت که از جمله کاربردهای ربات اسکارا نیز می‌باشد) روش‌های کلاسیک تنظیم ضرایب کنترل‌کننده‌های PD مانند زیگلر-نیکولز به خوبی جوابگو نیستند \cite{17}. در استفاده از روش‌هایی چون منطق فازی نیز، احتیاج به یک سری اطلاعات در مورد سیستم برای استخراج قوانین فازی می‌باشد \cite{18}. به منظور تنظیم بهینه ضرایب کنترل‌کننده‌های PD در این مقاله از روش یادگیری کیو افزایشی ارائه‌شده در \cite{19}، که به حداقل اطلاعات در مورد سیستم و محیط مورد تعامل نیاز دارد استفاده شده است.



\section{مدل دینامیکی ربات}

معادله کلی دینامیکی بازوی ربات $n$ لینکی را در نظر بگیرید.
\begin{equation}
M_q({q}) \ddot{{q}}+C_q({q},\dot{{q}})\dot{{q}}+g_q({q})=\tau-J^TF_{ext}
\end{equation} 

بردار $q\in \mathbb{R}^{(n_q \times 1)}$، ماتریس $M_q(q) \in \mathbb{R}^{(n_q \times n_q)}$ و $C_q(q) \in \mathbb{R}^{(n_q \times n_q)}$ به ترتیب نشان‌دهنده موقعیت مفاصل، ماتریس اینرسی و ماتریس کوریولیس ربات هستند. بردار $g_q({q})$ با ابعاد ${(n_q \times 1)}$ بیانگر اثر جاذبه و $\tau \in \mathbb{R}^{(n_q \times 1)}$ گشتاور کنترلی اعمال شده به ربات می‌باشد.  


\section{طراحی قانون کنترل امپدانس}

بلوک دیاگرام پیشنهادی الهام گرفته شده از \cite{20} (پایان نامه ملانی کیمل) در شکل (1) را در نظر بگیرید. هدف در حلقه کنترلی بالا تولید مسیر مطلوب از روی مسیر مرجع با توجه به نیروهای خارجی و موانع محیطی است.
\begin{figure*}[t]
 \centering
\includegraphics[height=8cm,width=13cm]{blockdiagram.png}
\caption{ بلوک دیاگرام پیشنهادی}
\label{TWENDYONE}
\end{figure*}
 در حلقه کنترلی پایین پس از تولید مسیر مطلوب در فضای مفصلی، با استفاده از کنترل‌کننده‌های PD که در آن ضرایب کنترل‌کننده‌ها از روش یادگیری کیو افزایشی تولید می‌شود، گشتاور کنترلی لازم محاسبه و به ربات اعمال می‌گردد.

 
  

فرض کنید ارتباط بین موقعیت مجری نهایی و متغیرهای مفصلی با تابع $f : \mathbb{R}^{(n_q \times 1)}\rightarrow \mathbb{R}^{(n_p \times 1)}$ بیان شود.
\begin{equation}
p(t)=f(q)
\end{equation}

با دو بار مشتق‌گیری از رابطه (2)، $\ddot{q}$ به صورت زیر محاسبه می‌شود.
\begin{equation}
\ddot{{q}}=J^{-1}(\ddot{p}-\dot{J}\dot{{q}})
\end{equation}



با جایگذاری $\ddot{q}$ در (1) داریم:

\begin{equation}
M_q({q}) J^{-1}(\ddot{p}-\dot{J}\dot{q})+C_q({q},\dot{{q}})\dot{{q}}+g_q({q})=\tau-J^TF_{ext}
\end{equation}
و با استفاده از قانون فیدبک خطی‌ساز و انتخاب $\tau$ به صورت زیر 

\begin{equation}\label{eq:5}
\begin{aligned}
\tau &= M_q({q}) J^{-1}(M_p^{-1}a_p-\dot{J}\dot{{q}})+C_q({q},\dot{{q}})\dot{{q}}  \\
& \qquad {} +g_q({q})+J^TF_{ext} 
\end{aligned}
\end{equation}

ترم‌های غیرخطی معادله (4) حذف و معادله به فرم خطی زیر در فضای کارتزین تبدیل می‌شود.
\begin{equation}
\ddot{p}=M_p^{-1}a_p
\label{eq11-1}
\end{equation}
اگر خطای ردیابی را به صورت زیر در نظر بگیریم:
\begin{equation}
e=p-p_{ref}
\end{equation}
با انتخاب $a_p$ به فرم 
\begin{equation}
\begin{aligned}
a_p &= M_p{\ddot{p}_{ref}}-D_p(\dot{p}-\dot{p}_{ref}) \\
& \qquad {} -K_p(p-p_{ref})+F_{ext} 
\end{aligned}
\end{equation}

امپدانس مطلوب بین خطای موقعیت و نیروهای خارجی به صورت زیر تعریف می‌گردد.
\begin{equation}
M_p\ddot{e}+D_p\dot{e}+K_pe=F_{ext}
\end{equation}
با جایگذاری $a_p$ انتخابی در رابطه (5) و ساده‌سازی، گشتاور کنترلی به صورت زیر محاسبه می‌گردد:
\begin{equation}
\begin{aligned}
 \tau &= \Lambda \ddot{p}_{ref}- \Lambda M_p^{-1}(D_p\dot{e}+K_pe-\dot{J}\dot{{q}})   \\
 & \qquad{} + (\Lambda M_p^{-1}+J^T)F_{ext} \\& \,+C_q({q},\dot{{q}})\dot{{q}}+g_q({q})
\end{aligned}
\end{equation}
که $\Lambda= M_q  J^{-1}$ است.

 \begin{figure*}[t]
 \centering
\includegraphics[height=5cm,width=13cm]{constraintregion2.png}
\caption{فضای قابل قبول حالت‌ها \cite{20}.}
\label{space}
\centering
\end{figure*}

 \subsection{انواع محدودیت‌ها}

\section{کنترل نامتغیر}

فرم کلی یک سیستم زمان پیوسته در رابطه (11) که $u$ ورودی سیستم می‌باشد را در نظر بگیرید.
 \begin{equation}
\dot{x}=f(x)+g(x)u 
\end{equation}

یک مجموعه $\mathcal{G}$، یک مجموعه مثبت نامتغیر\LTRfootnote{Positively Invariant Control} برای سیستم دینامیکی معادله (11) نامیده می‌شود اگر هر مسیر حالت سیستم که از یک نقطه در درون $\mathcal{G}$ آغاز می‌شود برای تمام زمان‌های بعدی در آن باقی بماند \cite{13}. 
\begin{equation*}
x(t) \in \mathcal{G}  \Rightarrow x(t+\Delta t) \in\mathcal{G} ,         \qquad  \forall   \Delta t> 0 
\end{equation*}   


با تغییر نمایش معادله دینامیکی (1) به صورت 
{
\begin{equation}
\begin{bmatrix}
\dot{q}\\ \ddot{q}
\end{bmatrix} 
=
\begin{bmatrix}
\dot{q} \\
-M_q^{-1}(C_q\dot{q}+g_q)
\end{bmatrix} 
+
\begin{bmatrix}
0 \\ 
M_q^{-1}
\end{bmatrix} 
(\tau -J^TF_{ext})
\end{equation}
}
آن را به فرم استاندارد معادله (11) در می‌آوریم. قانون $\tau$، ورودی کنترلی و $x=[q^T \,\, \dot{q}^T]^T$ حالت‌های سیستم هستند. در ادامه ابتدا نحوه مدل‌سازی انواع موانع محیطی و تعیین تابع محدودیت که به عنوان خروجی مجازی سیستم در نظر گرفته می‌شود، شرح داده خواهد شد. سپس با خطی‌سازی
 ورودی-خروجی\LTRfootnote{Input-Output Linearization} اثر سیگنال کنترلی بر نقض یا عدم نقض محدودیت بررسی و در نهایت با پیش‌بینی رفتار تابع محدودیت قانون کنترل نامتغیر طراحی می‌گردد. 
 
 محدودیت‌ها به طور کلی می‌توانند روی سیگنال کنترل‌، خروجی و یا متغیرهای حالت تعریف شوند. در این مقاله محدودیت‌ها روی متغیرهای حالت در نظر گرفته می‌شوند. برای برقراری ایمنی، ربات در محیط نباید از حد مشخصی به اجسام نزدیک‌تر شود. برای تعیین این حد مشخص ابتدا باید محدودیت‌ها مدل شوند. دو فرم معمول و کلی که برای مدل کردن اجسام در فضا به کار می‌روند، محدودیت‌های کروی\LTRfootnote{Spherical Constraints} و خطی هستند. یک دیوار در فضا را میتوان با محدودیت خطی مدل کرد. برای مدل‌سازی اجسام با شکل پیچیده‌تر میتوان آنها را به صورت ترکیبی از کره‌هایی که مکان مرکز کره‌ها متغیر با زمان است، در نظر گرفت.

تابع $h_c(x,\eta)=\eta_a^\mathsf{T}x-\eta_b$ بیانگر فاصله محدودیت خطی تا متغیرهای حالت یا همان مکان مجری نهایی ربات است. روابط (13) و (14) دو حالت خاص از این محدودیت که حد پایین یا بالای متغیرهای حالت را مدل می‌کنند، نشان می‌دهند.
 \begin{equation}
h_c(x,\eta)=\eta_{low}-x_i 
\end{equation}
\vspace{-1.2cm}
\begin{equation}
h_c(x,\eta)=x_i -\eta_{up}
\end{equation}
در مدل‌سازی محدودیت کروی با توجه به اینکه متغیرهای حالت داخل یا بیرون از کره قابل قبول هستند، می‌توان محدودیت را به فرم
\begin{equation}
h_c(x,\eta)=\eta_{r}-\norm{x-\eta_c}_{2}
\end{equation}
برای حالت‌های قابل قبول بیرون از کره و 
\begin{equation}
h_c(x,\eta)=\norm{x-\eta_c}_{2}-\eta_{r}
\end{equation}
برای داخل کره به کار برد. $\eta_c$ مکان مرکز کره و $\eta_r$ شعاع کره است. شکل (2) نواحی قابل قبول محدودیت‌های ذکر شده را نشان می‌دهد. مجموعه $l$  محدودیت متغیرهای حالت در یک بردار به صورت 
 
 \begin{equation}
y_c=h_c(x(t),\eta(t))=
\begin{bmatrix}
h_{c,1}(x(t),\eta_1(t)) \\
\vdots\\
h_{c,l}(x(t),\eta_l(t))
\end{bmatrix} 
\label{eq11-1}
\end{equation}
  نشان داده می‌شوند. تعداد محدودیت‌ها در مجموعه $\mathcal{B}$  ذخیره می‌شوند.
                                                                                                                                                                                                                                                                   \vspace{-.5cm}                                                                                                                                                                                                                                    

\begin{equation}
\mathcal{B}=\{ i \in \mathbb{N} \mid 1 \leq i \leq l \}
\end{equation}
  مجموعه قابل قبول متغیرهای حالت، متغیرهایی هستند که به ازای آنها هیچکدام از محدودیت‌ها نقض نشوند و یا به عبارت دیگر 
\vspace{-.3cm}
\begin{equation}
\mathcal{H}(t)=\{ x \in  \mathbb{R}^{2n_q}  \mid h_{c,i} \leq 0 \, \forall i \in \mathcal{B} \}
\end{equation}
برقرار باشد
 \subsection{خطی‌سازی ورودی-خروجی}
 
 
 
 به منظور طراحی قانون کنترل نامتغیر و جلوگیری از نقض محدودیت‌ها لازم است ابتدا اثر سیگنال کنترلی بر توابع محدودیت بررسی شود. هر محدودیت به عنوان یک خروجی مجازی برای سیستم به فرم استاندارد در معادله (11) در نظر گرفته می‌شود.
 \begin{equation} 
 \begin{cases} 
 \dot{x}=f(x)+g(x)u                   \\ 
  y_{c,i}=h_{c,i}(x(t),\eta(t))  
  \end{cases} 
\end{equation}
 
 خطی‌سازی ورودی-خروجی، سیستم در معادله (18) را با $2n_q$ حالت به یک زنجیره انتگرالی\LTRfootnote{Integral Chain} که خروجی آن همان $y_{c,i}$ اما ورودی $u$ با شبه ورودی $z_i$ جایگزین شده، منتقل می‌کند.
 \begin{equation} 
 z_i=y_{c,i}^{(r_i)}=a^\mathsf{T} _{c,i}(x,\eta)u + b_{c,i}(x,\eta,...,\eta^{(r_i)})
\end{equation}
 از آنجا که در سیستم‌های رباتیکی اغلب محدودیت‌ها روی سرعت مفاصل، مکان مفاصل و یا مجری نهایی اعمال می‌شوند، درجه نسبی سیستم $1$و یا $2$ خواهد بود. در بخش شبیه‌سازی محدودیت روی مکان مجری نهایی در نظر گرفته شده است و با دو بار مشتق‌گیری ورودی کنترلی $\tau$ ظاهر می‌شود و درجه نسبی سیستم $r_i=2$ خواهد بود. با دوبار مشتق‌گیری از $y_{c,i}(x(t),\eta(t))$ که حالت‌های $x(t)$ همان $q(t)$ها هستند و با توجه به معادله (12) داریم:
 
 \begin{align}
&\begin{aligned}
\dot{y}{_{c,i}} &= \dfrac{\partial h_{c,i}}{\partial {q}}\dot{q} +\dfrac{\partial h_{c,i}}{\partial \eta}\dot{\eta} 
\end{aligned}\nonumber\\ 
&\begin{aligned}
  \ddot{y}{_{c,i}}&=  \dfrac{\partial \dot{y}_{c,i}}{\partial \dot{q}}\ddot{q} + \dfrac{\partial \dot{y}_{c,i}}{\partial q}\dot{q}+ \dfrac{\partial \dot{y}_{c,i}}{\partial \eta}\dot{\eta} +\dfrac{\partial \dot{y}_{c,i}}{\partial \dot{\eta}}\ddot{\eta} \\
  & {} =\dfrac{{\partial h_{c,i}}}{{\partial q}}M^{-1}_q\tau -\dfrac{{\partial h_{c,i}}}{{\partial q}}M^{-1}_q(C_q\dot{q}+g_q)  \\
  & \qquad +\dfrac {\partial}{\partial q} \bigg(\dfrac{{\partial h_{c,i}}}{{\partial {q}}}\dot{q}\bigg)\dot{q}+2 \dfrac {\partial}{\partial \eta}\bigg(\dfrac{{\partial h_{c,i}}}{{\partial {q}}}\dot{q}\bigg)\dot{\eta} \\
  &\qquad +\dfrac {\partial}{\partial \eta}\bigg(\dfrac{{\partial h_{c,i}}}{{\partial {\eta}}}\dot{\eta}\bigg)\dot{\eta} +\dfrac{\partial h_{c,i}}{\partial \eta}\ddot{\eta} 
  \end{aligned}
\end{align}
 

بنابراین 
\begin{equation} 
{a}^\mathsf{T} _{c,i}=\dfrac{{\partial h_{c,i}}}{{\partial {q}}}M^{-1}_q
\end{equation}

\begin{align}
\begin{split}
b_{c,i} &= -\dfrac{{\partial h_{c,i}}}{{\partial q}}M^{-1}_q(C_q\dot{q}+g_q) +\dfrac {\partial}{\partial q} \bigg(\dfrac{{\partial h_{c,i}}}{{\partial {q}}}\dot{q}\bigg)\dot{q} \\
  &+2 \dfrac {\partial}{\partial \eta}\bigg(\dfrac{{\partial h_{c,i}}}{{\partial {q}}}\dot{q}\bigg)\dot{\eta}+\dfrac {\partial}{\partial \eta}\bigg(\dfrac{{\partial h_{c,i}}}{{\partial {\eta}}}\dot{\eta}\bigg)\dot{\eta} +\dfrac{\partial h_{c,i}}{\partial \eta}\ddot{\eta}
   \end{split}
\end{align}

است.

\subsection{پیش‌بینی رفتار تابع محدودیت}
هدف نهایی کنترل نامتغیر تصحیح سیگنال کنترل نامی در زمان‌هایی است که این سیگنال باعث نقض محدودیت‌ها می‌گردد. بنابراین قانون کنترل نامتغیر باید طوری انتخاب شود که حالت‌های سیستم را درون مجموعه قابل قبول حالت‌ها (مجموعه $\mathcal{H}$ ) نگه دارد. اگر شبه ورودی معادله (19) منفی انتخاب شود در نهایت باعث منفی شدن مقدار تابع محدودیت خواهد شد. به طور مثال اگر برای سیستم با درجه نسبی دو شبه ورودی یک عدد ثابت و منفی انتخاب شود یعنی مشتق دوم نیز منفی است. منفی بودن مشتق دوم باعث کاهش مشتق اول می‌شود. زمانی که مقدار مشتق اول منفی شود، مقدار تابع $y_{c,i}$ شروع به کاهش می‌کند و در نهایت به یک مقدار منفی می‌رسد. 

در صورتی که سیستم در لحظه $t$ باشد مقدار تابع $y_{c,i}$ بین لحظه $t$ و $t+\Delta t $ باید محاسبه شود. اگر مقدار ماکزیمم $y_{c,i}$ بین این دو بازه منفی باشد سیستم با سیگنال کنترل نامی به کار خود ادامه می‌دهد، اگر این مقدار صفر باشد سیگنال کنترل نامی باید اصلاح شود و سیستم با سیگنال اصلاحی به کار خود ادامه دهد و اگر مقدار ماکزیمم مثبت شود، به این معنی است که سیگنال اصلاحی قادر به نگه داشتن حالت‌های سیستم درون مجموعه قابل قبول حالت‌ها نیست.


با توجه به $\ddot{y}{_{c,i}(t)}=z_{i}$ در (19) و اینکه $z_{i}$ برابر با مقداری ثابت و منفی در نظر گرفته می‌شود داریم:
\vspace{-.5cm}

\begin{equation}
\dot{y}(t+\Delta t)=\dot{y}(t) +  \Delta t z_i
\end{equation}


از $\dot{y}$ به صورت زیر انتگرال می‌گیریم:
\vspace{-.4cm}

\begin{equation}
\int_{t}^{t+\Delta t}{\dot{y}_{c,i}(\tau)}d\tau = y_{c,i}(t+\Delta t)-y_{c,i}(t)
\end{equation} 
در نتیجه
\vspace{-.4cm}

\begin{equation}
y_{c,i}(t+\Delta t)=y_{c,i}(t)+ \int_{t}^{t+\Delta t}{\dot{y}_{c,i}(\tau)}d\tau 
\end{equation}
با تغییر متغیر $\tau=t+\chi$ و با توجه به معادله (25) خواهیم داشت:
\vspace{-.4cm}


\begin{equation}
y_{c,i}(t+\Delta t)=y_{c,i}(t)+ \int_{0}^{\Delta t}{\dot{y}_{c,i}(t+\chi)}d\chi 
\end{equation}

در نهایت با جایگذاری رابطه (23) در معادله (26) و ساده‌سازی، $y_{c,i}(t+\Delta t)$ به صورت زیر محاسبه می‌شود.
\begin{equation}
 y_{c,i}(t+\Delta t)  =y_{c,i} + \Delta t  \, \dot{y}_{c,i} +\dfrac{{(\Delta t)}^{2}}{2}z_i  
 \end{equation}
همانطور که بیان شد اگر شبه ورودی برای سیستم با درجه نسبی یک به صورت $z_i=\gamma_i\leq 0$ و برای سیستم با درجه نسبی بزرگتر از یک به صورت $z_i=\gamma_i< 0$ انتخاب شود باعث کاهش مقدار تابع خروجی $y_{c,i}$ می‌شود، بنابراین برای تابع خروجی یک مقدار ماکزیمم وجود دارد. برای تعیین لحظه $t$ که سوئیچ بین کنترل نامی و نامتغیر باید انجام گیرد، از تابع نامتغیر\LTRfootnote{Invariance Function} $\Phi_i(x,\eta,\gamma_i)$ استفاده می‌شود. تابع نامتغیر $\Phi_i(x,\eta,\gamma_i)$  مقدار ماکزیمم تابع محدودیت $y_{c,i}$  را بین زمان $t$ و $t+\Delta t$ نشان می‌دهد.

 برای سیستم با درجه نسبی 2، دو حالت مختلف پیش می‌آید. با مشتق از تابع خروجی $y_{c,i}(t+\Delta t)$ در (27) داریم:
 \begin{equation}
\dfrac{\partial y_{c,i}(t+\Delta t) }{\partial \Delta t}  =  \Delta t  \, \gamma_i   +  \dot{y}_{c,i}
 \end{equation}
اگر $ \dot{y}_{c,i}\leq 0$ باشد آنگاه عبارت (28) منفی است و مقدار فعلی تابع $y_{c,i}$، ماکزیمم است. در صورتی که $ \dot{y}_{c,i}>0$ باشد مقدار ماکزیمم در لحظه
  \begin{equation}
\Delta t =-\dfrac{\dot{y}_{c,i}}{\gamma_i  }  
 \end{equation}
     


اتفاق می‌افتد و بنابراین تابع نامتغیر برابر با مقدار زیر خواهد شد. 
\begin{equation}
\Phi_i(x,\eta,\dot{\eta},\gamma_i)=
 \begin{cases} 
  y_{c,i}                                                                      &                \dot{y}_{c,i} \leq 0  \\ 
  y_{c,i}-\dfrac{1}{2\gamma_i}\dot{y}{_{c,i}^2} &                \dot{y}_{c,i}>0
  \end{cases}
\end{equation}


\subsection{طراحی قانون کنترل نامتغیر}
 
در غیاب محدودیت‌ها سیستم با قانون کنترل نامی به کار خود ادامه می‌دهد. اگر سیگنال کنترلی را در این حالت $u_{no}$ بنامیم شبه ورودی برابر با مقدار زیر است.
\vspace{-.2cm}
\begin{equation} 
 z_{no,i}=a^\mathsf{T} _{c,i}(x,\eta)u_{no} + b_{c,i}\big(x,\eta,...,\eta^{(r_i)}\big)
\end{equation}
سوئیچ به کنترل اصلاحی باید زمانی انجام شود که حداقل یک محدودیت فعال باشد. محدودیت فعال به این معنی است که طبق پیش‌بینی یک از محدودیت‌ها قرار است نقض بشود. این محدودیت‌ها را با مجموعه $\mathcal{B}_{act}$ نشان می‌دهیم.
 \begin{equation}
\mathcal{B}_{act}=\{ i \in \mathcal{B}| \, \Phi_i(x,\eta,\dot{\eta},\gamma_i) \geq 0 \}
\end{equation}
 

مجموعه $\mathcal{G}$ به صورت
 \begin{equation}
\mathcal{G}=\{ x \in  \mathbb{R}^{2n_q}  | \Phi_i(x,\eta,\dot{\eta},\gamma_i) \leq 0 \quad  \forall i \in \mathcal{B} \}
\end{equation}
که مجموعه نامتغیر نامیده می‌شود شامل تمام متغیرهای حالتی است که به ازای آنها هیچکدام از محدودیت‌ها نقض نمی‌شوند و یا در آینده نیز منجر به نقض محدودیت‌ها نخواهند شد. زمانی که تابع $\Phi_i=0$ است بر روی مرز مجموعه نامتغیر قرار داریم.
\vspace{-.4cm}
 \begin{equation}
\partial \mathcal{G}=\{ x \in  \mathbb{R}^{2n_q}  | \Phi_i(x,\eta,\dot{\eta},\gamma_i) = 0 \,\,\,\, \forall i \in \mathcal{B} \}
\end{equation}

برای اینکه مجموعه $\mathcal{G}$ نامتغیر بماند یعنی هیچکدام از توابع $\Phi_i$ مثبت نشوند، باید زمانی که به مرز مجموعه رسیدیم تابع $\Phi_i$ شروع به کاهش کند. بنابراین مشتق تابع $\Phi_i$ به ازای متغیرهای حالت عضو $\partial \mathcal{G}$ باید منفی شود. مشتق تابع $\Phi_i$ بیان‌شده در (32) را به صورت زیر در نظر بگیرید:
\vspace{-.4cm}

 \begin{equation} 
 \dot{\Phi}_i(x,\eta,\dot{\eta},\gamma_i)=
 \begin{cases} 
    \dot{y}_{c,i}                                                                                &                \dot{y}_{c,i} \leq 0  \\ 
  \dot{y}_{c,i}(1-\dfrac{1}{\gamma_i}  \ddot{y}_{c,i})                      &                \dot{y}_{c,i}>0
  \end{cases} 
\end{equation}
 اگر $\dot{y}_{c,i} \leq 0$ باشد مقدار $\dot{\Phi}_i$ نیز منفی است و بنابراین تابع ${\Phi}_i$ کاهش پیدا می‌کند. اما  اگر $\dot{y}_{c,i}>0$ باشد در صورتی که $z_i=\ddot{y}_{c,i}>\gamma_i$  باشد تابع $\Phi_i$ در حال افزایش است. با توجه به مطالب بیان‌شده شبه ورودی تصحیح‌کننده برای سیستم با درجه نسبی دو به صورت زیر تعیین می‌گردد.
\vspace{-.5cm}

\begin{equation} 
 z_{c,i}=
 \begin{cases} 
    \gamma_i                  &            (z_{no,i} > \gamma_i    )  \wedge ({\Phi}_i \geq 0)  \\  
  z_{no,i}                       &            (z_{no,i} \leq \gamma_i )  \vee     ({\Phi}_i < 0) 
  \end{cases} 
\end{equation} 
با توجه به (21)، (37) و (38) برای اینکه مجموعه $\mathcal{G}$  نامتغیر بماند باید شرایط زیر برقرار باشد:
\begin{equation} 
 z_i \leq z_{c,i} \quad  \forall i \in \mathcal{B}_{act}
\end{equation}

  سیگنال کنترلی سیستم استاندارد معادله (20) باید طوری تعیین شود که در عین جلوگیری از نقض محدودیت‌های فعال، تاحد امکان به سیگنال کنترل نامی نزدیک باشد. با مینیمم‌سازی نرم اقلیدسی\LTRfootnote{Euclidean Norm} اختلاف سیگنال نامی و $u$ می‌توان به این هدف رسید \cite{14}.
\begin{equation}
\begin{aligned} 
u_c &= 
\argmin_{u} \norm{u-u_{no}}_{2}^2  \\ 
\text{s.t.} &\,\, a^\mathsf{T} _{c,i}(x,\eta)u + b_{c,i}(x,\eta,...,\eta^{(r_i)}) \leq z_{c,i} \,\,  \forall i \in \mathcal{B}_{act} \end{aligned}
\end{equation}    
  
  فرم کلی تابع هزینه را به صورت $f=\norm{Ax-B}_{2}^2$ در نظر بگیرید. مشتق اول $f=(Ax-B)^T(Ax-B)$ برابر با
  \begin{equation}
  \nabla f = 2A^T(Ax-B)
  \end{equation}
  و مشتق دوم به صورت زیر محاسبه می‌گردد.
  \begin{equation}
  \nabla ^2 f =2A^TA
  \end{equation}
  بنابراین ماتریس هسین\LTRfootnote{Hessian Matrix} تابع $f$ برابر با $H=2A^TA$ است. ماتریس هسین مثبت معین و بنابراین تابع $f$ اکیدا محدب می‌باشد. قیود مسأله در (40) نسبت به $u$ خطی و در نتیجه محدب هستند. با نوجه به محدب بودن $f$ و قیود، مسأله بهینه‌سازی (40) یک مسأله بهینه‌سازی محدب است و هر مینیمم محلی، مینیمم سراسری و یکتای مسأله می‌باشد \cite{21}.  
  
   مشتق تابع هزینه در مسأله (40) نسبت به $u$ به صورت $2u^T-2u^T_{no}$ را برابر صفر قرار می‌دهیم و بنابراین جواب مینیمم‌سازی زمانی که هیچ محدودیتی فعال نباشد به صورت زیر محاسبه می‌شود.
\vspace{-.5cm}
\begin{equation}
u=u_{no}
\end{equation}
چواب بهینه $u$ باید در محدودیت‌ها نیز صدق کند یعنی
\begin{equation}
u  \leq {A} ^{+}_{\mathcal{B}_{act}}(\bm{z}_{c,\mathcal{B}_{act}}-\bm{b}_{c,\mathcal{B}_{act}}) 
\end{equation}
که در
 آن ${\bm{z}}_{c, \mathcal{B}_{act}}=[z _{c,i}]_{i \in {\mathcal{B}_{act}}}$، ${A}_{\mathcal{B}_{act}}=[a^\mathsf{T} _{c,i}]_{i \in {\mathcal{B}_{act}}}$ 
 و ${\bm{b}}_{c, \mathcal{B}_{act}}=[a^\mathsf{T} _{c,i}]_{i \in {\mathcal{B}_{act}}}$ است، باید برقرار باشد. ماتریس شبه معکوس 
 مور-پنرز\LTRfootnote{Moor-Penroze} به ازای ${A}_{\mathcal{B}_{act}}$ به فرم زیر تعیین می‌گردد.
\begin{equation}
{A} ^{+}_{\mathcal{B}_{act}}={A} ^{T}_{\mathcal{B}_{act}}({A}_{\mathcal{B}_{act}}{A} ^{T}_{\mathcal{B}_{act}})^{-1}
\end{equation}   


با توجه اینکه در مسأله بهینه‌سازی محدب هر مینیمم محلی، مینیمم سراسری و یکتای مسأله است و با توجه به روابط (43) و (44) جواب تحلیلی مسأله بهینه‌سازی به صورت
\vspace{-.2cm}
\begin{equation}
u_c={A} ^{+}_{\mathcal{B}_{act}}(\bm{z}_{c, \mathcal{B}_{act}}-\bm{z}_{no, \mathcal{B}_{ act}})+ u_{no}
\end{equation} 
محاسبه می‌شود. ماتریس 
$\bm{z}_{no,\mathcal{B}_{ act}}$
 برابر با 
$[z_{no,i}]_{i \in {\mathcal{B}_{act}}}$
 است.

 
\section{طراحی کنترل‌کننده IQ-PD}


در این بخش ابتدا شرح مختصری از مفهوم یادگیری تقویتی\LTRfootnote{Reinforcement Learning}، اجزاء آن و یادگیری کیو داده خواهد شد. سپس نحوه محاسبه ضرایب کنترل‌کننده‌های PD در بلوک دیاگرام پیشنهادی شکل (1) با استفاده از یادگیری کیو افزایشی\LTRfootnote{Incremental Q-Learning} بیان می‌شود.

\subsection{یادگیری تقویتی}
%\linebreak
یادگیری تقویتی زیرمجموعه‌ای از یادگیری ماشین است که عامل یا همان یادگیرنده را قادر می‌سازد در محیط‌های پویا و تصادفی با استفاده از آزمون و خطا عملی را برای نزدیک‌تر شدن به هدف انجام دهد. عامل با استفاده از بازخورد اعمال انجام‌شده و تجربیات رفتار خود را بهبود می‌بخشد. یادگیری تقویتی با چهار عنصر اصلی سیاست\LTRfootnote{Policy}، تابع پاداش\LTRfootnote{Reward Function}، تابع ارزش\LTRfootnote{Value Function} و مدل محیط\LTRfootnote{Enviroment Model}، بر پایه‌ی فرآیندهای تصمیم‌گیری مارکوف\LTRfootnote{MDPs} بنا شده است.

توالی کار در فرآیندهای تصمیم‌گیری مارکوف به این صورت است که عامل و محیط در دنباله ای از گام‌های زمانی گسسته $t=0, 1, 2, 3, ...$  با هم در تعامل هستند. در هر گام 
زمانی $t$ عامل که در حالت ${x}_t$ قرار دارد عمل $k_t$ را انجام می‌دهد. بعد از انجام عمل، سیستم در زمان $t+1$ قرار گرفته و عدد اسکالر  $r_{w_{t+1}}$ را به عنوان پاداش عمل انجام شده دریافت می‌کند و سپس به حالت جدید $x_{t+1}$ منتقل می‌شود. مجموعه اطلاعات در MDPs به صورت $x_0, k_0, r_{w_1}, x_1, k_1, r_{w_2},... $ است.

هر فرآیند تصمیم‌گیری مارکوف با یک چهارتایی متشکل از فضای حالت‌ها $\mathbb{X}$، فضای عمل‌ها $\mathbb{K}$، احتمال انتقال حالت $\mathcal{P}$ و تابع پاداش $r_w$ مدل می‌شود. فرض کنید هر عملی که انجام می‌گیرد حالت $x_t$  به $x_{t+1}$ منتقل و سیگنال پاداش $r_{w_{t+1}}$ دریافت می‌شود. در صورتی که فرآیند تصمیم‌گیری یک فرآیند مارکوف باشد احتمال انتقال حالت $x_t$  به $x_{t+1}$ از $ \mathcal{P}\{x_{t+1}|x_t, x_{t-1}, ... , x_1, x_0\}= \mathcal{P}\{x_{t+1}|x_t\}$ پیروی می‌کند.
\begin{figure*}[t]
 \centering
\includegraphics[height=4cm,width=13cm]{updatingstate.jpg}
\caption{به روز رسانی فضای حالت در روش یادگیری کیو افزایشی \cite{20}.}
\label{space}
\centering
\end{figure*}

 
احتمال انتخاب عمل $k_t$ از طرف عامل زمانی که در حالت $x_t$ قرار دارد، با نگاشت $\pi:\mathbb{X} \to \mathbb{K} $  تعریف می‌شود و می‌توان نوشت:
\vspace{-.4cm}
\begin{equation}
 Pr\{k_t|x_t\}=\pi(x_t)
\end{equation}
نگاشت $\pi$ با نام سیاست شناخته می‌شود و هدف هر مسأله یادگیری تقویتی آن است که با انتخاب درست سیاست، پاداش تجمعی در طول زمان را حداکثر کند. برای مقایسه سیاست‌های مختلف می‌توان معیاری را برای سنجش آنها تعریف نمود. این معیار به عنوان خروجی سیاست در نظر گرفته می‌شود. خروجی یک سیاست، میزان پاداشی است که با اتخاذ تصمیمات متوالی با تبعیت از سیاست مذکور جمع آوری می‌شود. 
مقدار پاداش جمع‌آوری شده در طول زمان که در آن  $\gamma$ فاکتور تخفیف نامیده می‌شود
 برابر  
 $ R_t=r_{w_{t+1}} +\gamma \,r_{w_{t+2}} +\gamma^2 \,r_{w_{t+3}}+...=\sum_{t=0}^{n} \gamma^t \, r_{w_{t+1}}$ 
 است. هدف یادگیری، انتخاب سیاست بهینه است به طوری که پاداش تجمعی مورد نظر حداکثر شود. مفهوم مورد انتظار را با امید ریاضی نشان می‌دهیم. بنابراین تابع ارزش هر سیاست که خروجی آن سیاست است به فرم 
\begin{equation}
V^{\pi}(x_t)=E_{\pi}\{R_t|x_t\}=E_{\pi}\bigg \{\sum_{t=0}^{n} \gamma^t \, r_{w_{t+1}} \bigg \}
\end{equation}
بیان می‌گردد.
با ساده‌سازی مقدار پاداش جمع‌آوری شده و با توجه به رابطه (48) تابع ارزش به صورت زیر محاسبه می‌شود.
\vspace{-.3cm}
\begin{equation}
 V^{\pi}(x_t)=r_{w_{t+1}}+\gamma\, V^{\pi}(x_{t+1})
\end{equation}
سیاست بهینه که منجر به پاداش حداکثر می شود به صورت زیر خواهد بود.
\vspace{-.5cm}
\begin{equation}
V^{*}(x_t)= \argmax_{t}{V^{\pi}(x_t|k)} 
\end{equation}
  عمل $k$ همه عمل‌های ممکن به ازای $x_{t+1}$ است. بیان دیگر معادله (50) به صورت زیر است:
 \vspace{-.3cm}
\begin{equation}
 V^{*}(x_t)=r_{w_{t+1}}+\gamma\, V^{*}(x_{t+1}|k)
\end{equation}


با تغییر نمایش تابع ارزش، معادله (51) به (52) تبدیل می‌شود:
\vspace{-.5cm}
\begin{equation}
 Q^{*}(x_t,k_t)=r_{w_{t+1}}+\gamma\, Q^{*}(x_{t+1}, all\, actions)
\end{equation}


در مرجع \cite{22} ساده ترین فرم تابع ارزش یادگیری کیو به صورت زیر بیان شده است:
\begin{equation}
\begin{aligned}
Q(x_t,k_t) &= Q(x_t,k_t) + \alpha \big [ r_{w_{t+1}}  \\
& \quad{}+\gamma\, \max Q(x_{t+1}, \text{all actions})-Q(x_t,k_t) \big]  
\end{aligned}
\end{equation}
پارامتر $\alpha \in [0,1]$ نرخ یادگیری نامیده می‌شود.


\subsection{محاسبه ضرایب کنترل‌کننده‌های PD با استفاده از یادگیری کیو افزایشی}

کنترل‌کننده PD را در حالت پیوسته زمان در نظر بگیرید.
\begin{equation}
u(t)=k_pe(t) +k_d\dfrac{de(t)}{dt}
\end{equation}
خطای $e(t)$ برابر با اختلاف متغیرهای مفصلی و مقدار مطلوب آنها تعریف می‌شود.
\vspace{-.4cm}
\begin{equation}
e(t)=q-q_d
\end{equation}
با توجه به تعریف مشتق در فرم زمان گسسته به صورت $\dfrac{e(t_k)}{dt}=\dfrac{e(t_k)-e(t_{k-1})}{\Delta t}$  قانون کنترلی معادله (54) در فرم زمان گسسته با نرخ نمونه‌برداری $\Delta t$ به صورت زیر محاسبه می‌گردد.
\begin{equation}
u(t_k) = k_p e(t_k)+k_d \dfrac{e(t_k)-e(t_{k-1})}{\Delta t}
\end{equation}
  
  ضرایب کنترل‌کننده‌های PD و متغیرهای مفصلی ربات پیوسته و  این در حالی است که فضای حالت و عمل‌ها در یادگیری کیو افزایشی گسسته هستند. بنابراین ابتدا باید ضرایب PD و متغیرهای مفصلی ربات که به عنوان حالت و عمل انتخاب می‌شوند به نحو مناسبی گسسته‌سازی\LTRfootnote{Discretization} شوند. فضای حالت و عمل‌ها به صورت زیر نمایش داده می‌شوند.
\begin{equation}
\begin{cases}
{k}_t=(k_p^1,k^1_d, k^2_p,k^2_d, k^3_p, k^3_d, k^4_p, k^4_d) \\
x_t=(\theta_1,\theta_2,D_3,\theta_4)
 \end{cases}
\end{equation}  

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

\subsubsection{تولید فضای حالت}     
  
  
  شکل (3) که نحوه به روز رسانی فضای حالت را نشان می‌دهد در نظر بگیرید.  
  فضای حالت در زمان $t$ شامل سه حالت مشخص شده است. در زمان $t+1$ بررسی می‌کنیم که حالت جدید دریافت‌شده، $x'$،  در همسایگی کدام یک از حالت‌های موجود با شعاع $\rho$ است. چون در همسایگی هیچ حالتی نیست پس $x'$ به عنوان یک حالت جدید به فضای حالت اضافه می‌شود. اما در زمان $t+2$ چون حالت در همسایگی $x_2$ است فضای حالت تغییری نمی‌کند و به روز رسانی نمی‌شود. بنابراین در
لحظه $t+2$ فضای حالت برابر با $\mathbb{X}=\{x_0, x_1, x_2, x_3\}$ است. 
  
\subsubsection{حافظه موقت}       

حافظه موقت\LTRfootnote{Temporal Memory} که از چهارتایی $\mathcal{M}=(c_t, k_t, \ell_t, h_t)$ تشکیل شده است، وظیفه نگهداری موقت اطلاعات در حین فرآیند یادگیری را به عهده دارد. مرکز $c_t$ نام حالت‌های عضو $\mathbb{X}=\{c_1,c_2,...\}$ است. هر مرکز $c_t$ نماینده حالت‌های اطراف خود در فاصله اقلیدسی کمتر از $\rho$ می‌با‌شد. با توجه به (57) هر عمل $k_t$ متشکل از ضرایب کنترل‌های PD است. دیگر عنصر حافظه یعنی سطح گسسته‌سازی $0 \leq \ell_t \leq \mathcal{L}$ بیانگر سطح گسسته‌سازی فعلی در حین یادگیری می‌باشد و سطح کوانتیزاسیون $h_t$ مکان $i$ و $j$ مقدار ماکزیمم تابع ارزش به ازای یک حالت و همه عمل‌های ممکن برای آن نشان می‌دهد. 

فرض کنید $\xi$ یک متغیر باینری باشد. زمانی که حافظه بدون تغییر باقی می‌ماند مقدار آن برابر $1$ و زمانی که حافظه تغییر می‌کند $0$ خواهد بود. برای تعیین مقدار $\xi$، $N$ حافظه اخیر مقایسه می‌شوند. بنابراین اگر $\{\mathcal{M}_t=\mathcal{M}_{t-1}=...=\mathcal{M}_{t-N} , t \geq N\}$ برقرار باشد به این معنی است که حافظه بدون تغییر مانده و $\xi=1$ است.  

\subsection{گسسته‌سازی فضای عمل‌ها}
 \begin{figure*}[t]
 \centering
\includegraphics[height=7cm,width=15cm]{quanti2.png}
\caption{نحوه ایجاد سطح گسسته‌سازی و کوانتیزاسیون \cite{20}.}
\label{quantization}
\centering
\end{figure*}
در هر مرحله زمانی ربات (یادگیرنده) با استفاده از سیاست $\pi$، عمل‌ها یا همان ضرایب کنترل‌کننده‌های PD را انتخاب می‌کند. 
عمل $k_t=\{k^1,..., k^j,..., k^{\mathcal{J}}\}$ که $k^j$ برابر با $\{k_p^j,k_d^j\}$ است را در نظر بگیرید. اندیس $\mathcal{J}$  برابر با تعداد کنترل‌کننده‌های PD و بنابراین تعداد کل عمل‌ها $n_k=2\times\mathcal{J}$ است.

فرض کنید در حالت شروع، سیستم در سطح گسسته‌سازی $\ell_t=0$ قرار دارد. با ضرب کارتزین\LTRfootnote{Cartesian Product} تعداد حالت‌های موجود برای هر ضریب، تعداد $k^{\ell_t}$ عمل مختلف که برای ربات قابل انتخاب هستند، ایجاد می‌شود. در هر سطح گسسته‌سازی عمل‌ها با فاصله $\delta_{\ell_t} \in \Delta$ گسسته می‌شوند. برای روشن‌تر شدن نحوه گسسته‌سازی فرض کنید $\delta_0=\dfrac{1}{3}$ و ضرایب PD در شروع
 برابر $k_0=\bigg(\dfrac{1}{3},\dfrac{2}{3},\dfrac{3}{3},\dfrac{1}{3},\dfrac{1}{3}, \dfrac{1}{3}, \dfrac{1}{3},\dfrac{1}{3}\bigg) $  انتخاب شوند. اگر ضرایب بین 
 بازه $[0,1]$ اجازه تغییر داشته باشند چهار انتخاب برای هر ضریب وجود خواهد داشت. ضرب کارتزین آنها تعداد $4^8$ حالت مختلف برای ضرایب ایجاد می‌کند. ربات در هر سطح گسسته‌سازی و کوانتیزاسیون خاص می‌تواند عمل‌های ایجاد شده حول هر هسته گسسته‌سازی را انتخاب و امتحان کند. بنابراین در گسسته‌سازی فضای عمل‌ها مسأله اصلی، انتخاب هسته گسسته‌سازی است.
 
فرض کنید سیستم در زمان $t>N$ و سطح $\ell_t$ قرار دارد و عمل $k_t$ را با استفاده از روش $\epsilon-greedy$ از بین عمل‌های ممکن $k^{\ell_t}$ انتخاب کرده است. در هر گام زمانی تغییر حافظه‌ زمانی را چک می‌کنیم. اگر حافظه بدون تغییر مانده باشد ($\xi=1$) باید سطح گسسته‌سازی یکی افزایش یابد ($\ell_t'=\ell_t +1)$). در سطح گسسته‌سازی جدید عمل‌ها با فاصله $\delta_{l_t'}$ حول $k_t$ گسسته می‌شوند.

  هر بار که حافظه بدون تغییر بماند سطح گسسته‌سازی عوض شده و ممکن است یک سطح کوانتیزاسیون نیز ایجاد شود. اگر حالت و عملی را که به ازای آنها حافظه بدون تغییر مانده است $c_t$ و $k_t$ بنامیم زمانی که سطح گسسته‌سازی عوض می‌شود موقعیت $c_t$ و $k_t$ در جدول $Q$ به عنوان سطح کوانتیزاسیون ذخیره می‌شود.

  شکل (4) را در نظر بگیرید. برای سادگی و خواناتر شدن شکل در این بخش اندیس $t$ حذف شده است. در حین یادگیری، $N$ حافظه قبلی با مرکز و عمل $c_6$ و $k_4$ بدون تغییر مانده‌اند. بنابراین سطح گسسته‌سازی یکی افزایش داده می‌شود. سطح کوانتیزاسیون جدیدی با برچسب $[6,4]$ تشکیل می‌شود که هسته مرکزی عمل‌های آن $k_4$ است و برای ایجاد عمل‌های ممکن از این هسته با فاصله گسسته‌سازی $\delta_{\ell_t}$ استفاده می‌شود. در سطح گسسته‌سازی صفر فقط یک سطح کوانتیزاسیون وجود دارد. 
در هر سطح سطح گسسته‌سازی $\ell$ زیرفضاهایی\LTRfootnote{Subsections} از فضای حالت و عمل‌ها تشکیل می‌شود.

\columnbreak\noindent
\begin{minipage}{\columnwidth}
\centering
\setlength{\tabcolsep}{25pt}
\setLTRtable
\captionof{table}{پارامترهای مورد استفاده در شبیه‌سازی}
\begin{tabular}{cc}
\toprule[.4pt] \toprule[.4pt]
$m_1$ & $0.5$ \\ 
$m_2$ & $0.4$ \\ 
$m_3$ & $0.2$ \\ 
$m_4$ & $0.1$ \\ 
$I_{c_1}$ & $0.15$ \\ 
$I_{c_2}$ & $0.1$ \\ 
$I_{c_3}$ & $0.02$ \\ 
$I_{c_4}$ & $0.05$ \\ 
$M_p$ & $10*I _{3\times 3}(kg) $ \\ 
$D_p$ & $80*I _{3\times 3}(kg) $ \\ 
$K_p$ & r$600*I _{3\times 3}(kg) $ \\ 
\bottomrule[.4pt] \bottomrule[.4pt]
\end{tabular}
\end{minipage}

ادامه متن
\end{multicols}
\end{document}