کنترل بهینه گسسته

برق , مهندسی 36 بازدید

تعداد بازدید ها: 36

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

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

$$ \large X [ k + 1 ] = f ( x [ k ] , u [ k ] ) $$

که شرایط اولیه آن، $$ X[1] = X_0 $$ است.

می‌خواهیم کنترل $$ u[k] $$ را به گونه‌ای پیدا کنیم که تابع هزینه زیر را کمینه کند:

$$ \large J ( u ) = \Phi ( X _ f ) + \sum _ { k =1 } ^ { N } L ( X [ k ] , u [ k ] ) . $$

مشابه عملیاتی که برای یک مسئله پیوسته‌زمان در آموزش قبلی انجام دادیم، همیلتونین $$ H(X,u,\lambda ) $$ را به صورت زیر می‌نویسیم:

$$ \large H ( X [ k ] , u[ k ] , \lambda [ k + 1 ] ) = L ( X [ k ] , u [ k ] ) + \lambda [ k + 1 ] ^ T f ( x [ k ] , u [ k ] ) $$

در این حالت، شرایط بهینگی به صورت زیر خواهند بود:

$$ \large \left . \frac { \partial H } { \partial X } \right | _ { X [ k ] , u[ k ] , \lambda [ k + 1 ] } ^ T = \lambda [ k ] $$

$$ \large \left . \frac { \partial H } { \partial u } \right | _ { X [ k ] , u [ k ] , \lambda [ k + 1 ] } = 0 $$

$$ \large \left . \frac { \partial \Phi } { \partial X } \right | _ { X _ f } = \lambda _ f $$

رگولاتور خطی و گسسته درجه دوم با زمان محدود

معادله دینامیکی سیستم خطی زیر را در نظر بگیرید:

$$ \large X [ k + 1 ] = A x [ k ] + B u [ k ] $$

که شرایط اولیه آن $$ X[1] = X_0 $$ است.

هدف، یافتن $$ u[k] $$ است که تابع هزینه زیر را کمینه کند:

$$ \large J ( u ) = \frac { 1 } { 2 } X _ f ^ T S X _ f + \frac { 1 } { 2 } \sum _ { k = 1 } ^ { N } \left ( X [ k ] ^T Q X [ k ] + u [ k ] ^ T R u [ k ] \right ) . $$

با طی مراحلی مشابه آموزش کنترل بهینه، می‌توان نشان داد که شرایط لازم بهینگی به صورت زیر خواهند بود:

$$ \large \begin{align*}
u [ k ] & = – R ^ { – 1 } B ^ T \lambda [ k + 1 ] \\
\lambda [ k ] & = A ^ T \lambda [ k + 1 ] + Q X [ k ] \\
X [ k + 1 ] & = A x [ k ]+ B u [ k ] \end {align*} $$

و شرایط مرزی $$ X[1] = X_0 $$ و $$ \lambda[N+1] = S X_f $$ هستند.

از آنجا که معادلات در $$X$$ و $$u$$ خطی هستند، به دنبال جوابی به فرم زیر خواهیم بود:

$$ \large \lambda [ k ] = P [ k ] X [ k ] $$

با جایگذاری

$$ \large X[k+1] = Ax[k]+ B u[k] $$

$$\large u[k] = – R^{-1} B^T \lambda[k+1] = – R^{-1} B^T P[k+1] X[k+1] $$

در $$ u[k] = – R^{-1} B^T \lambda[k+1]$$، داریم:

$$ \large \begin{align*}
\lambda [ k + 1 ] & = P [ k + 1 ] X[ k + 1 ] = P [ k + 1 ] ( A X [k ] + B u [ k ] ) \\ & = P [ k + 1 ] ( A X [ k ] – B R ^ { – 1 } B ^ T \lambda [ k + 1 ] )
\end {align*} $$

$$ \large ( I + P [ k + 1 ] B R ^ { – 1 } B ^ T ) \lambda [ k + 1 ] = P [ k + 1 ] A X [ k ] $$

$$ \large \lambda [ k + 1 ] = ( I + P [ k + 1 ] B R ^ { – 1 } B ^ T ) ^ { – 1 } P [ k + 1 ] A X [ k ] $$

با جایگذاری در معادله انتشار کمک‌حالت، خواهیم داشت:

$$ \large \begin{align*}
P [ k ] X [ k ] & = A ^ T \left ( ( I + P [ k + 1 ] B R ^ { – 1 } B ^ T ) ^ { – 1 } P [ k + 1 ] A X [ k ] \right ) + Q X [ k ] \\
P [ k ] X [ k ] & = A ^ T ( I + P [ k + 1 ] B R ^ { – 1 } B ^ T ) ^ { – 1 } P[ k + 1 ] A X [ k ] + Q X [ k ] \end {align*} $$

از آنجایی که معادله بالا برای همه $$ X[k]$$ها برقرار است، داریم:

$$ \large P [ k ] = A ^ T ( I + P [ k + 1 ] B R ^ { – 1 } B ^ T ) ^ { – 1 } P [ k +1 ] A + Q $$

فرم دیگر

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

$$ \large u [ k ] = – R ^ { – 1 } B ^ T \lambda [ k + 1 ] $$

$$ \large R u [ k ] = – B ^ T P [ k + 1 ] X [ K + 1 ] $$

$$ \large R u [ k ] = – B ^ T P [ k + 1 ] ( A X [ k ] + B u [ k ] ) $$

$$ \large ( R + B ^ T P [ k + 1 ] B ) u [ k ] = – B ^ T P [ k + 1 ] A X [ k ] $$

$$ \large u [ k ] = – ( R + B ^ T P [ k + 1 ] B ) ^ { – 1 } B ^ T P [ k + 1] A X [ k ] $$

معادلات کمک‌حالت نیز برابرند با:

$$ \large \lambda [ k ] = A ^ T \lambda [ k + 1 ] + Q X [ k ] $$

$$ \large \lambda [ k ] = A ^ T P [ k + 1 ] X [ k + 1 ] + Q X [ k ] $$

$$ \large \lambda [ k ] = A ^ T P [ k + 1 ] ( A X [ k ] + B u [ k ] ) + Q X [ k ] $$

$$ \large P [ k ] X [ k ] = A ^ T P [ k +1 ] ( A – B ( R + B ^ T P [ k + 1 ] B ) ^ { – 1 } B ^ T P [ k + 1 ] A ) X [ k ] + Q X [ k ] $$

از آنجایی که رابطه بالا برای همه $$X[k] $$ها برقرار است، داریم:

$$ \large P [ k ] = A ^ T P [ k + 1 ] A – A ^ T P [ k + 1 ] B ( R + B ^ T P [ k + 1 ] B ) ^ { – 1 } B ^ T P [ k + 1 ] A + Q $$

رگولاتور خطی و گسسته درجه دوم با زمان نامحدود

معادله یک سیستم خطی دینامیکی را با شرایط اولیه $$X[1] = X_0 $$ به صورت زیر در نظر بگیرید:

$$ \large X [ k + 1 ] = A x [ k ]+ B u [ k ] $$

هدف، یافتن $$ u[k] $$ است که تابع زیر را کمینه کند:

$$ \large J ( u ) = \sum _ { k = 1 } ^ { \infty } \left ( X [ k ] ^ T Q X [ k ] + u [ k ] ^ T R u [ k ] \right ) . $$

تحت شرایط حالت ماندگار یا فرض زمان بی‌نهایت، $$ P[k] = P[k+1] $$ و معادله ریکاتی به رابطه زیر کاهش می‌یابند:

$$ \large P = A ^ T + A ^ T ( I + P B R ^ { – 1 } B ^ T ) ^ { – 1 } P A + Q $$

با توجه به فرمول دوم، معادله ریکاتی را می‌توان به صورت زیر نوشت:

$$ \large P = A ^ T P A – A ^ T P B ( R + B ^ T P B ) ^ { – 1 } B ^ T P A + Q $$

توجه کنید که دو فرم معادله ریکاتی معادل هستند و با استفاده از لم معکوس ماتریس می‌توان هر یک را به دیگری تبدیل کرد.

جواب معادله ریکلاتی جبری گسسته را می‌توان با استفاده از دستور “dare” در متلب پیاده‌سازی کرد:

clc
close all
clear all 

t_f = 20;
dt = 0.0001;
t_dis = 0:dt:t_f;
N = length(t_dis);
S{N} = 100*eye(2);
R = dt;
Q = .0001*eye(2)*dt;
A = [1 dt;0 1];
B = [0 ;dt];
K{N} = [0 0];

for i = N-1:-1:1
    %K{i} = inv(R + B'*S{i+1}*B)*B'*S{i+1}*A;
    %S{i} = Q + K{i}'*R*K{i} + (A-B*K{i})'*S{i+1}*(A-B*K{i}); % First form
    %S{i} = Q+ A'*S{i+1}*A - A'*S{i+1}*B*inv(R+B'*R*B)*B'*S{i+1}*A;% Second form    
    S{i} =  A'*inv(eye(2) + S{i+1}*B*inv(R)*B')*S{i+1}*A + Q; % Third form
    K_norm(i) = norm(K{i});
end

X(:,1) = [1;0];
X_dlqr(:,1) = [1;0];
P_dlqr = dare(A,B,Q,R);
K_dlqr = inv(R)*B'*P_dlqr;

for i = 2:N
    u(i-1) = -inv(R)*B'*S{i-1}*X(:,i-1);
    X(:,i) = A * X(:,i-1) + B*u(i-1);
    X_dlqr(:,i) = A * X_dlqr(:,i-1) - B*K_dlqr*X_dlqr(:,i-1) ;

end


figure;
subplot(2,1,1)
plot(t_dis,X(1,:),t_dis,X_dlqr(1,:))
legend('DyP','LQR')
ylabel('position')
subplot(2,1,2)
plot(t_dis,X(2,:),t_dis,X_dlqr(2,:))
legend('DyP','LQR')
ylabel('velocity')
xlabel('time')

clc

close all

clear all

t_f = 20;

dt = 0.0001;

t_dis = 0:dt:t_f;

N = length(t_dis);

S{N} = 100*eye(2);

R = dt;

Q = .0001*eye(2)*dt;

A = [1 dt;0 1];

B = [0 ;dt];

K{N} = [0 0];

for i = N-1:-1:1

%K{i} = inv(R + B'*S{i+1}*B)*B'*S{i+1}*A;

%S{i} = Q + K{i}'*R*K{i} + (A-B*K{i})'*S{i+1}*(A-B*K{i}); % First form

%S{i} = Q+ A'*S{i+1}*A - A'*S{i+1}*B*inv(R+B'*R*B)*B'*S{i+1}*A;% Second form

S{i} = A'*inv(eye(2) + S{i+1}*B*inv(R)*B')*S{i+1}*A + Q; % Third form

K_norm(i) = norm(K{i});

end

X(:,1) = [1;0];

X_dlqr(:,1) = [1;0];

P_dlqr = dare(A,B,Q,R);

K_dlqr = inv(R)*B'*P_dlqr;

for i = 2:N

u(i-1) = -inv(R)*B'*S{i-1}*X(:,i-1);

X(:,i) = A * X(:,i-1) + B*u(i-1);

X_dlqr(:,i) = A * X_dlqr(:,i-1) - B*K_dlqr*X_dlqr(:,i-1) ;

end

figure;

subplot(2,1,1)

plot(t_dis,X(1,:),t_dis,X_dlqr(1,:))

legend('DyP','LQR')

ylabel('position')

subplot(2,1,2)

plot(t_dis,X(2,:),t_dis,X_dlqr(2,:))

legend('DyP','LQR')

ylabel('velocity')

xlabel('time')

نتیجه اجرای برنامه بالا، نمودارهای زیر خواهد بود:

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

تذکر: جواب LQR با استفاده از توابع “care” و “dare” فقط به مسائلی با زمان نامحدود قابل اعمال هستند.

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

همه توابع هزینه را نمی‌توان به صورت مسئله LQR بیان کرد. حالتی را در نظر بگیرید که می‌خواهیم زمان را کمینه کنیم.
فرمول‌بندی‌های LQR برای شرایطی که قیود اضافه به سیستم اعمال شود به خوبی تعمیم‌یافته نیستند.
- برای قیود نامساوی کنترل، جواب کنترل LQR به مقادیری خاصی محدود می‌شود.
- کار با حالت یا ترکیب حالت-کنترل دشوارتر است و شرایط بهینگی منتجه بسیار پیچیده خواهد بود.

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

به عنوان حامی، استارتاپ، محصول و خدمات خود را در انتهای مطالب مرتبط مجله فرادرس معرفی کنید.

به اشتراک بگذارید:

منبع Vivek Yadav

سید سراج حمیدی

«سید سراج حمیدی» دانش‌آموخته مهندسی برق است. فعالیت‌های کاری و پژوهشی او در زمینه سیستم‌های فتوولتائیک و کاربردهای کنترل در قدرت بوده و، در حال حاضر، آموزش‌های مهندسی برق و ریاضیات مجله فرادرس را می‌نویسد.

بر اساس رای 1 نفر

آیا این مطلب برای شما مفید بود؟

مهندس در وکا (WEKA) ابزاری برای یادگیری ماشین و داده کاوی — راهنمای شروع به کار نوشته است:

متشکر
مجید عوض زاده در اگزرژی (Exergy) در مکانیک — به همراه مثال نوشته است:

با سلام و تشکر بابت بازخورد ارائه شده؛ برای طراحی بهینه یک سیستم، بهتر است آنتروپی تولیدی کل سیستم میبینمم شود؛ در این صورت اگزرژی دریافتی، ماکزیمم بوده و می‌توان به نقطه بهینه طراحی دست یافت.
علیرضا در اگزرژی (Exergy) در مکانیک — به همراه مثال نوشته است:

مقالتون عالی بود فقط سوالی که داشتم چطور میتونم یک سیستم مکانیکی رو از نظر اگزرژی بهینه کنم

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