حداقل مربعات در متلب

حداقل مربعات در متلب را با چند مثال مختلف بررسی می کنیم.

برای آموزش شبیه سازی روش های حداقل مربعات و حداقل مربعات بازگشتی و… در متلب مثال های عملی را با هم بررسی می کنیم.

فرض کنیم چنین سوالاتی مطرح است:

1-سیستم زیر را در نظر بگیرید که در آن ) e(tنویز سفید با واریانس 0/1است.

فرض کنید پارامتر a=0.7و b=2باشد با اعمال روش هاي
الف- حداقل مربعات و ب- حداقل مربعات بازگشتی تخمین پارامترهاي aو bرا بدست آورید و رسم نمایید.تغییرات پارامترها
را نسبت به تعداد تکرارها رسم نموده و نحوه همگرایی تخمین پارامترها را توضیح دهید.سپس روشهاي مختلف را با هم
مقایسه کنید.
-2فرض کنید مدل ریاضی سیستم بصورت زیر است

ورودي uبه سیستم و خروجی yبه صورت زیر است پارامترهاي a0و b0و b1را با روش حداقل مربعات تخمین بزنید.

-3سیستم

با ورودي PRBSکه داراي دامنه ±1است در نظر بگیرید با فرض ) e(tنویز سفید با واریانس واحد شبیه سازي کنیددو
مدل براي سیستم زیر در نظر بگیرید

با در نظر گرفتن N=100این مدلها را با روش هاي
حداقل مربعات
حداقل مربعات تعمیم داده شده
متغییرهاي کمکی
تخمین زده و با هم مقایسه کنید

با پروژه آماده متلب که در زیر قرار داده شده و آموزش آن با کامنت گذاری مشخص شده است بخوبی میتوان همه مراحل را دید:

=================================

clc;پاک کردن صفحه متلب

clear all;پاک کردن متغیرهای گذشته و بلا استفاده

close all;بستن نمودارهای باز

%==========================================================================

%Q1سوال شماره 1

a=0.7;مقدار پارامتر

b=2;مقدار پارامتر

sigm_e=0.1;واریانس نویز

N=100;تعداد داده ها

N1=N;متغیر کمکی برای ذخیره تعداد داده ها

y=zeros(1,N);تعریف بردار خروجی

u=50*idinput(N);ورودی شناسایی

ep=wgn(1,N,0);نویز سفید

sigm_ep=1/(N)*sum(ep.^2);واریانس نویز

e=sqrt(sigm_e/(sigm_ep))*ep;ایجاد نویز با واریانس مورد نظر

جمع آوری داده از سیستم

for t=2:N

y(t)=-a*y(t-1)+b*u(t-1)+e(t);

end

y1=y;u1=u;

ترسیم داده های شناسایی

figure (1)

subplot(2,1,1);stairs(1:1:N,y,’b’,’linewidth’,2);grid on;axis([1 N -400 400]);ylabel(‘y’);

subplot(2,1,2);stairs(1:1:N,u,’b’,’linewidth’,2);grid on;axis([1 N -60 60]);ylabel(‘u’);xlabel(‘sample’);

%Least squaresروش حداقل مربعات

Phi=zeros(N-1,2);ماتریس داده ها

Y=zeros(N-1,1);بردار خروجی

for t=2:N

Phi(t,:)=[-y(t-1),u(t-1)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;محاسبه پارامترها

teta1=teta;

نمایش مقادیر داده در صفحه متلب

disp(‘Q1: Least squares estimation for N=100’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b=’);disp(teta(2));

N=200;تکرار شبیه سازی برای تعداد داده های بیشتر

N2=N;

y=zeros(1,N);

u=50*idinput(N);

ep=wgn(1,N,0);

sigm_ep=1/(N)*sum(ep.^2);

e=sqrt(sigm_e/(sigm_ep))*ep;

for t=2:N

y(t)=-a*y(t-1)+b*u(t-1)+e(t);

end

y2=y;u2=u;

Phi=zeros(N-1,2);

Y=zeros(N-1,1);

for t=2:N

Phi(t,:)=[-y(t-1),u(t-1)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;

teta2=teta;

disp(‘Q1: Least squares estimation for N=200’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b=’);disp(teta(2));

N=300;تکرار شبیه سازی برای تعداد داده های بیشتر

N3=N;

y=zeros(1,N);

u=50*idinput(N);

ep=wgn(1,N,0);

sigm_ep=1/(N)*sum(ep.^2);

e=sqrt(sigm_e/(sigm_ep))*ep;

for t=2:N

y(t)=-a*y(t-1)+b*u(t-1)+e(t);

end

y3=y;u3=u;

Phi=zeros(N-1,2);

Y=zeros(N-1,1);

for t=2:N

Phi(t,:)=[-y(t-1),u(t-1)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;

teta3=teta;

disp(‘Q1: Least squares estimation for N=300’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b=’);disp(teta(2));

N=400;تکرار شبیه سازی برای تعداد داده های بیشتر

N4=N;

y=zeros(1,N);

u=50*idinput(N);

ep=wgn(1,N,0);

sigm_ep=1/(N)*sum(ep.^2);

e=sqrt(sigm_e/(sigm_ep))*ep;

for t=2:N

y(t)=-a*y(t-1)+b*u(t-1)+e(t);

end

y4=y;u4=u;

Phi=zeros(N-1,2);

Y=zeros(N-1,1);

for t=2:N

Phi(t,:)=[-y(t-1),u(t-1)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;

teta4=teta;

disp(‘Q1: Least squares estimation for N=400’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b=’);disp(teta(2));

ترسیم درصد خطای نسبی تخمین پارامترها بر حسب تعداد داده ها

Nv=[100,200,300,400];

Er_Nv_a=[abs(a-teta1(1))/a,abs(a-teta2(1))/a,abs(a-teta3(1))/a,abs(a-teta4(1))/a]*100;

Er_Nv_b=[abs(b-teta1(2))/b,abs(b-teta2(2))/b,abs(b-teta3(2))/b,abs(b-teta4(2))/b]*100;

figure (2)

subplot(2,1,1);bar(Nv,Er_Nv_a,’b’);xlabel(‘N’);ylabel(‘a relative error (%)’);

subplot(2,1,2);bar(Nv,Er_Nv_b,’b’);xlabel(‘N’);ylabel(‘b relative error (%)’);

%Recursive least squaresروش حداقل مربعات بازگشتی

lmbd=0.98;ضریب فراموشی

P=1e6*eye(2); ماتریس

teta_rls1=zeros(2,N1);ماتریس تخمین ها

تخمین حداقل مربعات بازگشتی

for t=2:N1

phi=[-y1(t-1);u1(t-1)];

Kt=P*phi/(lmbd*eye(1)+phi.’*P*phi);

P=(eye(2)-Kt*phi.’)*P/lmbd;

teta_rls1(:,t)=teta_rls1(:,t-1)+Kt*(y1(t)-phi.’*teta_rls1(:,t-1));

end

ترسیم پارامترهای تخمینی به روش حداقل مربعات و مقادیر واقعی و مقادیر روش قبل

figure (3)

subplot(2,1,1);stairs(1:N1,teta1(1)*ones(1,N1),’b’,’linewidth’,2);grid on;

hold on;

stairs(1:N1,teta_rls1(1,:),’g’,’linewidth’,2);

hold on;

stairs(1:N1,a*ones(1,N1),’r’,’linewidth’,2);

legend(‘LS for N=100′,’RLS for N=100′,’a=0.7’);

subplot(2,1,2);stairs(1:N1,teta1(2)*ones(1,N1),’b’,’linewidth’,2);grid on;

hold on;

stairs(1:N1,teta_rls1(2,:),’g’,’linewidth’,2);

hold on;

stairs(1:N1,b*ones(1,N1),’r’,’linewidth’,2);

legend(‘LS for N=100′,’RLS for N=100′,’b=2’);

%==========================================================================

%Q2سوال شماره 2

u=[1,0.8,0.6,0.4,0.2,0,0.2,0.4,0.6,0.8,1,0.8,0.6,0.4,0.2];

y=[0.9,2.5,2.4,1.3,1.2,0.8,0,0.9,1.4,1.9,2.3,2.4,2.3,1.3,1.2];

N=length(u);

ترسیم داده های شناسایی

figure (4)

subplot(2,1,1);stairs(1:1:N,y,’b’,’linewidth’,2);grid on;axis([1 N -1 4]);ylabel(‘y’);

subplot(2,1,2);stairs(1:1:N,u,’b’,’linewidth’,2);grid on;axis([1 N -1 2]);ylabel(‘u’);xlabel(‘sample’);

Phi=zeros(N-1,3);

Y=zeros(N-1,1);

for t=2:N

Phi(t,:)=[-y(t-1),u(t),u(t-1)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;تخمین پارامترها

disp(‘Q2: Least squares estimation’);

disp(‘estimated a1=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

%==========================================================================

%Q3سوال شماره 3

پارامترها

a=-0.9;

b0=1;

b1=0.5;

c=1;

sigm_e=1;

N=200;

y=zeros(1,N);

u=idinput(N);

ep=wgn(1,N,0);

sigm_ep=1/(N)*sum(ep.^2);

e=sqrt(sigm_e/(sigm_ep))*ep;

for t=3:N

y(t)=-a*y(t-1)+b0*u(t-1)+b1*u(t-2)+c*e(t);

end

figure (5)

subplot(2,1,1);stairs(1:1:N,y,’b’,’linewidth’,2);grid on;axis([1 N -20 20]);ylabel(‘y’);

subplot(2,1,2);stairs(1:1:N,u,’b’,’linewidth’,2);grid on;axis([1 N -2 2]);ylabel(‘u’);xlabel(‘sample’);

% Model A — Least squares تخمین جداقل مربعات برای مدل اول

Phi=zeros(N-1,3);

Y=zeros(N-1,1);

for t=3:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;

disp(‘Q3: Least squares estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

%Model A — Extended Least squares تخمین حداقل مربعات تعمیم یافته برای مدل اول

Phi=zeros(N-1,4);

Y=zeros(N-1,1);

teta0=[teta(1);teta(2);teta(3);1];

teta=teta0;

for t=3:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2),y(t-1)-Phi(t-1,:)*teta];

Y(t)=y(t);

if (t<10)

teta=teta0;

else

teta=(Phi(1:t,:).’*Phi(1:t,:))\Phi(1:t,:).’*Y(1:t);

end

end

disp(‘Q3: Extended Least squares estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

disp(‘estimated c=’);disp(teta(4));

%Model A — Instrumental variables روش متغیرهای کمکی

Phi=zeros(N-1,3);

Z=zeros(N-1,3);

Y=zeros(N-1,1);

for t=4:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2)];

Z(t,:)=[u(t-1),u(t-2),u(t-3)]; متغیرهای کمکی انتخاب شده

Y(t)=y(t);

end

teta=(Z.’*Phi)\Z.’*Y;تخمین پارامترها

disp(‘Q3: Instrumental variable estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

% Model B — Least squares روش حداقل مربعات برای مدل دوم

Phi=zeros(N-1,3);

Y=zeros(N-1,1);

for t=3:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2)];

Y(t)=y(t);

end

teta=(Phi.’*Phi)\Phi.’*Y;

disp(‘Q3: Least squares estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

%Model B — Extended Least squares روش حداقل مربعات تعمیم یافته برای مدل دوم

Phi=zeros(N-1,5);

Y=zeros(N-1,1);

teta0=[teta(1);teta(2);teta(3);1;1];

teta=teta0;

for t=3:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2),y(t-1)-Phi(t-1,:)*teta,y(t-2)-Phi(t-2,:)*teta];

Y(t)=y(t);

if (t<10)

teta=teta0;

else

teta=(Phi(1:t,:).’*Phi(1:t,:))\Phi(1:t,:).’*Y(1:t);

end

end

disp(‘Q3: Extended Least squares estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

disp(‘estimated c0=’);disp(teta(4));

disp(‘estimated c1=’);disp(teta(5));

%Model B — Instrumental variables روش متغیرهای کمکی برای مدل دوم

Phi=zeros(N-1,5);

Z=zeros(N-1,5);

Y=zeros(N-1,1);

teta0=[teta(1);teta(2);teta(3);1;1];

teta=teta0;

for t=6:N

Phi(t,:)=[-y(t-1),u(t-1),u(t-2),y(t-1)-Phi(t-1,:)*teta,y(t-2)-Phi(t-2,:)*teta];

Z(t,:)=[u(t-1),u(t-2),u(t-3),y(t-4),y(t-5)]; متغیرهای کمکی

Y(t)=y(t);

if (t<10)

teta=teta0;

else

teta=(Z(1:t,:).’*Phi(1:t,:))\Z(1:t,:).’*Y(1:t);

end

end

disp(‘Q3: Instrumental variables estimation’);

disp(‘estimated a=’);disp(teta(1));

disp(‘estimated b0=’);disp(teta(2));

disp(‘estimated b1=’);disp(teta(3));

disp(‘estimated c0=’);disp(teta(4));

disp(‘estimated c1=’);disp(teta(5));

نگران نباشید ما این کد آماده متلب را در انتهای همین پست برای شما قرار داده ایم.

نتایج شبیه سازی را در زیر قرار داده ایم:

برای دانلود کد آماده این آموزش روی گوگل +1 کلیک کنید

شاید علاقه مند باشید مطالب دیگر سایت متلبی را نیز ببینید:

• پروژه کامپیوتری درس فرایند تصادفی با متلب
• طراحی و شبیه سازی درایو موتور رلوکتانس سوییچی با متلب
• شبیه سازی رفتار یک رله مغناطیسی با متلب
• سفارش شبیه سازی مقالات درس کنترل توان راکتیو
• تشخیص فونم ها(لب خوانی) با SVM در متلب
• شناسایی سیستم غیر خطی ربات بازوی مسطح دو درجه آزادی توسط شبکه عصبی
• شبیه سازی سیستم درایو کرامر استاتیکی با متلب
• حل تابع با سری تیلور و روش نیوتن در متلب
• توابع نمایی در متلب
• مشکل لایسنس متلب در زمان اجرای نرم افزار متلب
• پیش بینی سری زمانی به کمک شبکه عصبی در متلب