حالا اینها را که به دست آوردی بذار داخل یه بردار به نام a و ضرب در 2:10 کن
a.*(2:10)
تا مخرج به دست بیاد
صورت را هم که برنامه اش را گذاشتم
حالا اینها را که به دست آوردی بذار داخل یه بردار به نام a و ضرب در 2:10 کن
a.*(2:10)
تا مخرج به دست بیاد
صورت را هم که برنامه اش را گذاشتم
اینطوری شد :
clc; close all; clear;
a=[1.934 ,1.554 ,1.174 ,0.794 , 0.414 ,0.077012 ,0.057655 ,0.044775 ,0.035774];
soorat=a.*(2:10)
makhraj = zeros(1,10)
for numplayers=2:10
lb=zeros(1,numplayers);
ub=ones(1,numplayers);
[x1, fx1]=fmincon(@fffu,ones(1,numplayers)*.5,[],[],[],[],lb,ub);
disp('makhraj=')
disp(-fx1)
makhraj(numplayers) =(-fx1)
efficiency(numplayers) = soorat(numplayers-1)/(-fx1)
end
plot(2:10,efficiency(2:10).*(2:10)), grid
ylabel('Symmetric NE Efficiency');
xlabel('Number of Suppliers in Microgrid')
نمودار :
http://www.uplooder.net/img/image/83...2f7ae/fig3.png
این الآن درسته؟
آخرین ویرایش به وسیله coronaa : چهارشنبه 26 شهریور 1393 در 12:20 عصر
تو این برنامه که سایز صورت و مخرج یکی نیست!
مهندس وقتی ما payoff تغذیه کننده های مختلف رو بدست آوردیم دیگه نیازی نیست که ما بیایم مثلن برای تغذیه کننده های با تعداد 5 تا رو ضرب در 5 بکنیم.(این خودش payoff پنج تا تغذیه کننده هست!) منظورم این خطه :
soorat=a.*(2:10)
من الآن فقط خود بردار payoff رو بعنوان صورت در نظر گرفتم نمودارش اینطوری شد :
http://www.uplooder.net/img/image/84...c6eeb/fig5.png
معقول تره به نظرم!
اگه شکل (2) مقاله رو نگاه کنید دقیقن مثل نمودارهای upperbound و simulated شد.
آخرین ویرایش به وسیله coronaa : چهارشنبه 26 شهریور 1393 در 13:16 عصر
علتش اینه که توی برنامه ای که گذاشتی دو بار ضربدر 2:10 کردی
a.*2:10
و
plot(2:10,efficiency(2:10).*(2:10)), grid
که لازمه موقع plot دیگه در 2:10 ضرب نشه
پس برنامه کلی به این صورت شد با نمودار دومی که گذاشتم؟
clc; close all; clear;
payoff=[1.934 ,1.554 ,1.174 ,0.794 , 0.414 ,0.077012 ,0.057655 ,0.044775 ,0.035774];
soorat=payoff.*(2:10)
% soorat=payoff;
makhraj = zeros(1,10) ;
efficiency = zeros(1,10);
for numplayers=2:10
lb=zeros(1,numplayers);
ub=ones(1,numplayers);
[x1, fx1]=fmincon(@fffu,ones(1,numplayers)*.5,[],[],[],[],lb,ub);
disp('makhraj=')
disp(-fx1)
makhraj(numplayers) =(-fx1);
efficiency(numplayers) = soorat(numplayers-1)/makhraj(numplayers);
end
خب ما تا حالا از جدول 3 و 4 استفاده نکردیم. اینطوری باشه که این واسه همه ی مسائل که یکی در میاد!
سایز صورت و مخرج یکینیست. چیکار کنم؟
تا اینجا که پیش رفتیم مربوط به رقابت بین تولید کننده ها بود از این به بعد که در جدول 3و 4 اومده مربوط به همکاری بین اونهاست که فکر کنم از shapley استفاده کرده ( فکر کنم مقاله را به صورت کامل مطالعه نکردی)
مطالعه کردم ولی فکر میکردم که صورت کسر باید با نقطه تعادل نش بدست بیاد و مخرج کسر با گیم تئوری مثل shapley !
الآن باید چیکار کنیم؟
این قسمت رو که دیگه توی همون مقاله اولی بدست آوردیم دیگه. داریم کدش رو. با همون mat-tug
ما توی این برنامه تونستیم تازه مطابق شکل(4)، نمودار simulated efficiency رو رسم کنیم. نمودار upperbound و lowerbound هم از همین برنامه باید بدست بیاد؟
قبلن توی اولین کد یه همچین نموداری بدست آوردی که شبیه lowerbound هست. همینه یا اون کد اشتباس؟
http://www.uplooder.net/img/image/93...768b7/fig5.png
مهندس نمودار upperbound و lowerbound از روی همین برنامه بدست نمیآد؟
این کد زیر رو به برنامه اضافه کردم. نمودار relative error برای upperbound بدست اومد. شکل(4)
error(numplayers) = abs(ub(numplayers-1)-efficiency)/efficiency;
http://www.uplooder.net/img/image/16...upperbound.png
ولی lowerbound بدست نمیآد. فقط یه خط راست میشه. شاید ورودی lowerbound که همه رو صفر در نظر گرفتیم اشتباه باشه!
http://www.uplooder.net/img/image/60...d_relative.png
شما نظرتون چیه؟
آخرین ویرایش به وسیله coronaa : پنج شنبه 27 شهریور 1393 در 15:52 عصر
upperbound و lowerbound در ص 425 و در Corollary 2 توضیح داده شده
ببینم چی کار می کنی
من الآن روی مقدار lowerbound مشکل دارم. فکر میکنم صفر نباید باشه!
نمودار error of upperbound درسته؟ تائید کردی؟
یه مشکل دیگه ای هم که دارم. توی برنامه اصلی سایز صورت و مخرج یکی نیست!
من از این برنامه استفاده کردم. شما اون برنامه اولی رو گفتی اشتباهه!
clc; close all; clear;
payoff=[1.934 ,1.554 ,1.174 ,0.794 , 0.414 ,0.077012 ,0.057655 ,0.044775 ,0.035774];
soorat=payoff.*(2:10)
makhraj = zeros(1,10) ;
efficiency = zeros(1,10);
error =zeros(1,10);
for numplayers=2:10
lb=zeros(1,numplayers);
ub=ones(1,numplayers);
[x1, fx1]=fmincon(@fffu,ones(1,numplayers)*.5,[],[],[],[],lb,ub);
disp('makhraj=')
disp(-fx1)
makhraj(numplayers) =(-fx1);
efficiency(numplayers) = soorat(numplayers-1)/makhraj(numplayers);
error(numplayers) = abs(ub(numplayers-1)-efficiency)/efficiency;
end
plot(2:10,efficiency(2:10)), grid
ylabel('Symmetric NE Efficiency');
xlabel('Number of Suppliers in Microgrid')
axis([2 10 0 1.1])
% ----------------------------------------------
figure
plot(2:10,error(2:10)), grid
ylabel('relative error of the upperbound');
xlabel('Number of Suppliers in Microgrid')
% legend('cos_x','sin_x')
axis([2 10 0 1.1])
این lb,ub برابر با upperbound و lowerbound نیست که
خواهشا و لطفا به همون صفحه مراجعه کن
ای بابا ینی این مقادیر payoff و چیزهای دیگه رو بدست آوردیم همه اشتباه بود؟
پس نمودارا چرا درست در اومد؟؟؟
چشم! فقط منو نزن :))
من یه لحظه داشتم ناامید می شدم. فکر کردم همه اشتباهه.
باشه من میخونمش. ببینم میتونم کاری انجام بدم یا نه.
ممنونم.
سلام دوست عزیز
من یه فانکشن ساختم به اسم Lmax و یکی هم به اسم Jmax ، طبق روابط 16 و 17 . بصورت زیر :
function ret = Lmax(L)
global k;global D;global alpha;
ret=-((alpha-k)./m).*L.^2+k.*D.*L;
end
function ret = Jmax(J)
global k;global D;global alpha;
ret=-(((alpha./m-k).*J.^2)+K.*D.*J);
end
بعد اومدم توی برنامه اصلی حد بالا و پایین راندمان رو بصورت شرطی قرار دادم.(بالای روابط 16)
حالا دیگه نمیدونم مقادیر fminbnd و fmincon رو چطوری مشخص کنم.
لاندا از کجا باید بدست بیاد؟
clc; close all; clear;
global q;global k;global D;global alpha ;global lambda;
alpha=.3; beta=1.2; k=.38; D=6.3; m=10;
L=sum(q);
J=sum(q);
for numsupplier=2:10
if lambda > k*D/((m+1)*k-2*alpha)
lobnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/(((alpha-k).*Lmax^2/m)+k*D*Lmax);
upbnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/((alpha/m-k)*Jmax^2+k*m*D*Jmax);
end
if lambda < k*D/((m+1)*k-2*alpha)
lobnd= (m^2.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k).*Lmax.^2)+k*m*d*Lmax);
upbnd= (m.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k*m).*Jmax.^2)+k*m*d*Jmax);
end
[x1, fx1]=fmincon(@fff,Q,[],[],[],[],lobnd,upbnd);
for tekrar=1:400
for supplier=1:numsupplier
[x2, fx2]=fminbnd(@fff,0,1);
Q(supplier)=x2;
end
end
efficiency(numsupplier)=fx2/fx1;
end
plot(2:10,efficiency(2:10).*(2:10)), grid
ylabel('Symmetric NE efficiency');
xlabel('Number of Suppliers in Microgrid')
در ضمن، ظاهرن با روابط 40 و 41 و 42 میشه بازده رو بدست آورد. هم رابطه ی V_NE مشخصه و هم رابطه ی V_opt
اینجا دیگه انگار payoff هم نمیخواد!
چه دلیلی داشته که از این حلقه ها استفاده کردی؟
for tekrar=1:400
for supplier=1:numsupplier
بخشهایی را خوب پیش رفتی
اولا توضیح بده اون صفحه چیکار می خواهد بکنه و مراحلش را بگو
ثانیا در هلپ متلب در مورد دو تابع fmincon و fminbnd مطالعه کن و بگو هر کدوم برای چه کاری مناسبه و اینکه میشه یکی را به جای دیگری استفاده کرد؟
ما میخوایم حد بالا و پایین راندمان رو بدست بیاریم. یعنی مشخص کنیم که راندمانی که بدست میاد از این مقادیر نباید تجاوز کنه.
توی این کرانهای بالا و پایین، دو متغیر Lmax و Jmax هم دارم که باید مشخص بشن. یعنی خود این دوتا هم باید اول ماکسیمم بشن تا توی کرانها قرار بگیرن.
L و J مجموع مقادیر انرژی تغذیه کننده ها هستن که ما اول باید این دوتا رو با روابط 16 و 17 بیشینه کنیم و بعد توی کرانهای بالا و پایین قرار بدیم.
ضمن اینکه خود این راندمای که میخواد بدست بیاد بصورت شرطی باید باشه. یعنی اگه لاندا از اون مقدار بزرگتر بود ، کران بالا و پایین راندمان به این صورته و اگر لاندا کوچکتر از اون مقدار بود ، راندمان باید توی این بازه قرار بگیره. یعنی حد بالا و پایین یا همون upperbound و lowerbound آن به این صورت میشه.
ـــــــــــــــــــــــــ ـــــــــــــــــــــــــ ـــــــــــــــــــــــــ ــــــــــــــــــــــــ
در مورد توابع fmincon و fminbnd ، خدمت شما عرض کنم که هر دو برای مسائل غیر خطی مقید شده استفاده میشن. ولی fminbnd برای یک متغیر استفاده شده و fmincon برای توابع چند متغیره استفاده میشه. پس نمیشه جای همدیگه هم استفاده کرد!
به نظرم توی این مسئله هم باید از fmincon استفاده کنیم.
شاید بشه برای مسائل یک متغیره هم از fmincon استفاده کرد فقط به جای چند متغیر از یک متغیر استفاده می کنیم ( من که متلب ندارم شما امتحان کن ببین میشه یا نه)
سوال: در اینجا یک تابع داریم که در اون تابع باید همزمان دو تا متغیر max بشن
یا اینکه دو تا تابع داریم که برای هرکدوم باید یک متغیر max بشه؟
ببینم چی کار می کنی نمودار را بکش
طبق روابط 16 و 17 دو تا تابع داریم که باید ماکزیمم بشن و هر تابع هم یک متغیر جدا داره. L و J دوتا متغیر هستن که L برای کران پایین و J برای کران بالا استفاده میشه.
مهندس من سوادم در همین حده! خب همون برنامه ای که بالا گذاشتم در همین حد بلدم. ران که میکنم یه سری متغیرها مثل لاندا رو نمیشناسه.
باید چیکار کنم؟
آخرین ویرایش به وسیله coronaa : جمعه 28 شهریور 1393 در 09:54 صبح
یه همچین چیزی
بعضی ها را کد خودت کپی کردم ( ضمانتش با خودت)
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
for m=2:10
parameterfunL=@(L,apha,k,m,D) (alpha-k)/m*(L^2)+k*D*L;
parameterfunJ=@(j,apha,k,m,D) (alpha/m-k)*(j^2)+k*D*j;
fl=@(L) parameterfunL(L,apha,k,m,D);
fj=@(j) parameterfunJ(j,apha,k,m,D);
Lmax=fminbnd(@fl,0,m*lambda);
Jmax=fminbnd(@fj,0,m*lambda);
if lambda > k*D/((m+1)*k-2*alpha)
lobnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/(((alpha-k).*Lmax^2/m)+k*D*Lmax);
upbnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/((alpha/m-k)*Jmax^2+k*m*D*Jmax);
end
if lambda < k*D/((m+1)*k-2*alpha)
lobnd= (m^2.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k).*Lmax.^2)+k*m*d*Lmax);
upbnd= (m.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k*m).*Jmax.^2)+k*m*d*Jmax);
end
end
برای خط 9 و 10 اشاره گر رو برداشتم. خطا میداد :
Lmax=fminbnd(fl,0,m*lambda);
Jmax=fminbnd(fj,0,m*lambda);
خب حالا این upbnd و lobnd رو بدست میاده. یک عدد هست. همین درسته یا باید یک بردار بده؟
بالایی اشتباه شدم. به ازای m های مختلف up,lo بدست میاد ولی نمیدونم چطوری باید توی یه بردار بذارم!
حل شد.
الآن به دوتا بردار تبدیل شد ولی عدداش یه طوریه!
upbnd
1.0e+004 *
Columns 1 through 7
0 1.6698 1.5303 0.7856 0.7421 0.2988 0.0941
Columns 8 through 10
0.0000 0.0000 0.0000
lobnd
1.0e+004 *
Columns 1 through 7
0 3.3397 4.5908 3.1424 3.7103 1.7930 0.6588
Columns 8 through 10
0.4396 0.3148 0.2507
مهندس دوتا برنامه رو با هم ترکیب کردم و بجای lp و ub از برنامه جدید استفاده کردم تا ببینم نمودارش چطوری میشه.
این کارم اشتباهه؟
فانکشن fffu
function ret = fffu(q)
alpha=.3; beta=1.2; k=.38; D=6.3;
ret=-sum(alpha.*q.^2+k.*(D-sum(q)).*q);
end
برنامه اصلی :
clc; close all; clear;
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
for m=2:10
parameterfunL=@(L,alpha,k,m,D) (alpha-k)/m*(L^2)+k*D*L;
parameterfunJ=@(j,alpha,k,m,D) (alpha/m-k)*(j^2)+k*D*j;
fl=@(L) parameterfunL(L,alpha,k,m,D);
fj=@(j) parameterfunJ(j,alpha,k,m,D);
Lmax=fminbnd(fl,0,m*lambda);
Jmax=fminbnd(fj,0,m*lambda);
if lambda > k*D/((m+1)*k-2*alpha)
lobnd(m)=((m*(k^2)*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/(((alpha-k).*Lmax^2/m)+k*D*Lmax);
upbnd(m)=((m*(k^2)*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/((alpha/m-k)*Jmax^2+k*m*D*Jmax);
end
if lambda < k*D/((m+1)*k-2*alpha)
lobnd(m)= (m^2.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k).*Lmax.^2)+k*m*D*Lmax);
upbnd(m)= (m.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k*m).*Jmax.^2)+k*m*D*Jmax);
end
end
fprintf('Lowerbound\n ');
disp(lobnd)
fprintf('\n---------------------------------------\n');
fprintf('Upperbound \n');
disp(upbnd)
% ---------------------------------------------------
% ---------------------------------------------------
payoff=[1.934 ,1.554 ,1.174 ,0.794 , 0.414 ,0.077012 ,0.057655 ,0.044775,0.035774];
soorat=payoff.*(2:10);
makhraj = zeros(1,10) ;
efficiency = zeros(1,10);
error =zeros(1,10);
for numplayers=2:10
lb=lobnd(1:numplayers);
ub=upbnd(1:numplayers);
[x1, fx1]=fmincon(@fffu,ones(1,numplayers)*.5,[],[],[],[],lb,ub);
disp('makhraj=')
disp(-fx1)
makhraj(numplayers) =(-fx1);
efficiency(numplayers) = soorat(numplayers-1)/makhraj(numplayers);
error(numplayers) = abs(ub(numplayers-1)-efficiency)/efficiency;
end
plot(2:10,efficiency(2:10)), grid
ylabel('Symmetric NE Efficiency');
xlabel('Number of Suppliers in Microgrid')
axis([2 10 0 1.1])
% ----------------------------------------------
figure
plot(2:10,error(2:10)), grid
ylabel('relative error of the upperbound');
xlabel('Number of Suppliers in Microgrid')
% legend('cos_x','sin_x')
axis([2 10 0 1.1])
% ----------------------------------------------
یه پیغام خطا هم میده :
Exiting due to infeasibility: 1 lower bound exceeds the corresponding upper bound.
makhraj=
??? In an assignment A(I) = B, the number of elements in B and
I must be the same.
Error in ==> plot5 at 42
makhraj(numplayers) =(-fx1);
آخرین ویرایش به وسیله coronaa : جمعه 28 شهریور 1393 در 12:04 عصر
ظاهرن توی مرحله یک ، fx1 پوچ در میاد!
چیکار باید بکنم؟
آخرین ویرایش به وسیله coronaa : جمعه 28 شهریور 1393 در 16:06 عصر
مهندس من وقتی برای خط پایینی مقدار lb و ub رو بصورت زیر انتخاب میکنم، واسه اولین مرحله( یعنی m برابر 2)، fx1 رو پوج میده. این کارم اشتباهه؟
lb=lobnd(1:numplayers);
ub=upbnd(1:numplayers);
[x1, fx1]=fmincon(@fffu,ones(1,numplayers)*.5,[],[],[],[],lb,ub);
میخواستم طبق شکل (4)، نمودار حد بالا و پایین رو بدست بیارم.
خواهش میکنم یه کمکی بکن.
ین یکی را امتحان کن
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
ll=zeros(1,10);
uu=zeros(1,10);
for m=2:10
parameterfunL=@(L,alpha,k,m,D) -(alpha-k)/m*(L^2)+k*D*L;
parameterfunJ=@(j,alpha,k,m,D) -(alpha/m-k)*(j^2)+k*D*j;
fl=@(L) parameterfunL(L,alpha,k,m,D);
fj=@(j) parameterfunJ(j,alpha,k,m,D);
Lmax = fmincon(@fl,m*lambda/2,[],[],[],[],0,m*lambda);
Jmax = fmincon(@fj,m*lambda/2,[],[],[],[],0,m*lambda);
if lambda > k*D/((m+1)*k-2*alpha)
lobnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/(((alpha-k).*Lmax^2/m)+k*D*Lmax);
upbnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/((alpha/m-k)*Jmax^2+k*m*D*Jmax);
end
if lambda < k*D/((m+1)*k-2*alpha)
lobnd= (m^2.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k).*Lmax.^2)+k*m*D*Lmax);
upbnd= (m.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k*m).*Jmax.^2)+k*m*D*Jmax);
end
ll(m)=lobnd;
uu(m)=upbnd;
end
من توی خط 12 و 13 اشاره گرها رو بر میدارم. وگرنه خطا میده:
جواب :
ll =
0 Inf Inf Inf Inf Inf Inf Inf Inf Inf
>> uu
uu =
1.0e+015 *
Columns 1 through 7
0 Inf Inf 2.2085 Inf Inf Inf
Columns 8 through 10
0.0542 Inf Inf
دو تا فایل برای این تابعها درست کن
function ret=fj(j)
global m;
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
ret= -(alpha/m-k)*(j^2)+k*D*j;
end
function ret=fl(L)
global m;
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
ret= -(alpha-k)/m*(L^2)+k*D*L;
end
بعدش هم این یکی را اجرا کنچ
global m;
alpha=.3; beta=1.2; k=.38; D=6.3;lambda=1;
ll=zeros(1,10);
uu=zeros(1,10);
for m=2:10
Lmax = fmincon(@fl,m*lambda/2,[],[],[],[],0,m*lambda);
Jmax = fmincon(@fj,m*lambda/2,[],[],[],[],0,m*lambda);
if lambda > k*D/((m+1)*k-2*alpha)
lobnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/(((alpha-k).*Lmax^2/m)+k*D*Lmax);
upbnd=((m*k^2*D^2*(k-alpha))/(((m+1)*k-2*alpha))^2)/((alpha/m-k)*Jmax^2+k*m*D*Jmax);
end
if lambda < k*D/((m+1)*k-2*alpha)
lobnd= (m^2.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k).*Lmax.^2)+k*m*D*Lmax);
upbnd= (m.*(alpha.*lambda.^2+k.*(D-m*lambda).*lambda))/(((alpha-k*m).*Jmax.^2)+k*m*D*Jmax);
end
ll(m)=lobnd;
uu(m)=upbnd;
end
جواب :
http://uplood.ir/53lR