differential evalution code Error using * Inner matrix dimensions must agree.

1 view (last 30 days)
Error in @(x)sum((minus((x(1)),V)*sin(2.*pi.*x(2).*t+x(3)).^2))/n Error in noise_de (line 56) pop(i).Cost=CostFunction(pop(i).Position);
  1. clear all; close all; clc
  2. fs=200; %sampling freq.
  3. dt =1/fs;
  4. n=fs/3 %number of samples/cycle
  5. m=3 %no. of cycles
  6. fi=50;
  7. t = dt*(0:400); %data window
  8. ww=wgn(201,1,-40);
  9. size(transpose(ww))
  10. t =dt*(0:200);
  11. y=sin(2*pi*fi*t + 0.3);
  12. for j=0:200/(n*m)
  13. t =dt*(j*m*n:(j+1)*m*n);
  14. x=sin(2*pi*fi*t + 0.3)+transpose(wgn(1+n*m,1,-40));
  15. V=x
  16. tmax=0.01;
  17. lastreported=0;
  18. %% Problem Definition
  19. t_est=[];
  20. f_est=[];
  21. dt=1/fs;
  22. i_max=tmax*fs
  23. for ii=0:i_max
  24. if(ii/i_max*100-lastreported>=1)
  25. lastreported=ii/i_max*100;
  26. fprintf('%5.2f%%\n',lastreported);
  27. end
  28. t=(ii:ii+n-1)*dt;
  29. CostFunction=@(x) sum((minus((x(1)),V)*sin(2*pi.*x(2).*t+x(3)).^2))/n; % Cost Function
  30. nVar=3; % Number of Decision Variables
  31. VarSize=[1 nVar]; % Decision Variables Matrix Size
  32. VarMin=[0,48,0]; % Lower Bound of Decision Variables
  33. VarMax=[1000,52,2*pi]; % Upper Bound of Decision Variables
  34. %% DE Parameters
  35. MaxIt=200; % Maximum Number of Iterations
  36. nPop=50; % Population Size
  37. beta=0.5; % Scaling Factor
  38. pCR=0.2; % Crossover Probability
  39. minCost=1e-10;
  40. %% Initialization
  41. empty_individual.Position=[];
  42. empty_individual.Cost=[];
  43. BestSol.Cost=inf;
  44. pop=repmat(empty_individual,nPop,1);
  45. for i=1:nPop
  46. pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
  47. pop(i).Cost=CostFunction(pop(i).Position);
  48. if pop(i).Cost<BestSol.Cost
  49. BestSol=pop(i);
  50. end
  51. end
  52. BestCost=zeros(MaxIt,1);
  53. %% DE Main Loop
  54. for it=1:MaxIt
  55. for i=1:nPop
  56. x=pop(i).Position;
  57. A=randperm(nPop);
  58. A(A==i)=[];
  59. a=A(1);
  60. b=A(2);
  61. c=A(3);
  62. % Mutation
  63. %beta=unifrnd(beta_min,beta_max);
  64. y=pop(a).Position+beta.*(pop(b).Position-pop(c).Position);
  65. y = max(y, VarMin);
  66. y = min(y, VarMax);
  67. % Crossover
  68. z=zeros(size(x));
  69. j0=randi([1 numel(x)]);
  70. for j=1:numel(x)
  71. if j==j0 || rand<=pCR
  72. z(j)=y(j);
  73. else
  74. z(j)=x(j);
  75. end
  76. end
  77. NewSol.Position=z;
  78. NewSol.Cost=CostFunction(NewSol.Position);
  79. if NewSol.Cost<pop(i).Cost
  80. pop(i)=NewSol;
  81. if pop(i).Cost<BestSol.Cost
  82. BestSol=pop(i);
  83. end
  84. end
  85. end
  86. % Update Best Cost
  87. BestCost(it)=BestSol.Cost;
  88. % Show Iteration Information
  89. %disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
  90. if(minCost>BestSol.Cost)
  91. break;
  92. ErrorTarget=0.00000001;
  93. EvalMax=10000*n;
  94. end
  95. end
  96. %% Show Results
  97. % disp(['Iteration ' num2str(ii) ': Best Cost = ' num2str(BestSol.Position(2))]);
  98. t_est=[t_est;(ii)*dt];
  99. f_est=[f_est;BestSol.Position(2)];
  100. if(minCost>BestSol.Cost)
  101. %break;
  102. ErrorTarget=0.00000001;
  103. EvalMax=10000*n;
  104. end
  105. end
  106. end
  107. t_est
  108. f_est
  109. plot (t_est,f_est,'red')
  110. hold on
  111. xlabel('time')
  112. ylabel('frequency')
  113. title('DE white noise ')
  114. c=vpa(rms(fi(t_est)-f_est))
  115. plot (t_est,fi*ones(size(t_est)))
  116. hold off

Accepted Answer

Stephan
Stephan on 2 Dec 2020
@(x)sum((minus((x(1)),V).*sin(2.*pi.*x(2).*t+x(3)).^2))/n
% ^
% |
% ------- Elementwise multiplication
  19 Comments
Stephan
Stephan on 3 Dec 2020
Edited: Stephan on 3 Dec 2020
The problem is here:
t_est=[t_est;(ii+n-1)*dt];
this leads to values of t_est ranging from 0...1 - not 0...200 as expected. If you change it to:
t_est=[t_est; (ii-1)*dt*fs];
The values are scaled 0...200. But due to missing insight in the problem i dont know if this is correct. At least the plot seems to be ok.
Since fi is ascalar, it wont plot - i suggest to use yline instead:
plot(t_est,f_est,'red')
xlim([0 200])
yline(fi)
xlabel('time')
ylabel('frequency')
title('DE white noise')
So once again the whole code:
clear all; close all; clc
fs=200; %sampling freq.
dt =1/fs;
n=fs/3 %number of samples/cycle
m=3 %no. of cycles
fi=50;
t = linspace(0,200,1+n*m); %data window
v=@(t) (sin(2*pi*fi.*t + 0.3)+ transpose(wgn(1+n*m,1,-40)));
tmax=1;
t_est=[];
f_est=[];
dt=1/fs;
i_max=tmax*fs;
for ii=0:i_max
V=v(t);
CostFunction=@(x) sum((V-x(1)*sin(2*pi*x(2)*t+x(3))).^2)/n; % Cost Function
nVar=3; % Number of Decision Variables
VarSize=[1 nVar]; % Decision Variables Matrix Size
VarMin=[0,48,0]; % Lower Bound of Decision Variables
VarMax=[1000,52,2*pi]; % Upper Bound of Decision Variables
%% DE Parameters
MaxIt=1000; % Maximum Number of Iterations
nPop=50; % Population Size
beta=0.5; % Scaling Factor
pCR=0.2; % Crossover Probability
minCost=1e-10;
%% Initialization
empty_individual.Position=[];
empty_individual.Cost=[];
BestSol.Cost=inf;
pop=repmat(empty_individual,nPop,1);
for i=1:nPop
pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
pop(i).Cost=CostFunction(pop(i).Position);
if pop(i).Cost<BestSol.Cost
BestSol=pop(i);
end
end
BestCost=zeros(MaxIt,1);
%% DE Main Loop
for it=1:MaxIt
for i=1:nPop
x=pop(i).Position;
A=randperm(nPop);
A(A==i)=[];
a=A(1);
b=A(2);
c=A(3);
% Mutation
%beta=unifrnd(beta_min,beta_max);
y=pop(a).Position+beta.*(pop(b).Position-pop(c).Position);
y = max(y, VarMin);
y = min(y, VarMax);
% Crossover
z=zeros(size(x));
j0=randi([1 numel(x)]);
for j=1:numel(x)
if j==j0 || rand<=pCR
z(j)=y(j);
else
z(j)=x(j);
end
end
NewSol.Position=z;
NewSol.Cost=CostFunction(NewSol.Position);
if NewSol.Cost<pop(i).Cost
pop(i)=NewSol;
if pop(i).Cost<BestSol.Cost
BestSol=pop(i);
end
end
end
% Update Best Cost
BestCost(it)=BestSol.Cost;
% Show Iteration Information
%disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
if(minCost>BestSol.Cost)
break;
end
end
%% Show Results
disp(['Iteration ' num2str(ii) ': Best Cost = ' num2str(BestSol.Position(2))]);
t_est=[t_est; (ii-1)*dt*fs];
f_est=[f_est;BestSol.Position(2)];
end
t_est
f_est
RMSE = sqrt(mean((f_est-fi).^2))
plot(t_est,f_est,'red')
xlim([0 200])
yline(fi)
xlabel('time')
ylabel('frequency')
title('DE white noise')

Sign in to comment.

More Answers (0)

Categories

Find more on Graph and Network Algorithms in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!