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)

Tags

Community Treasure Hunt

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

Start Hunting!