How to tune this model with this PSO coding? Is there any mistake in this coding or do I need to add something?

1 view (last 30 days)
Below is the NARX integrated with PSO but i cant tune PSO with NARX.
clear;
%load input-output data
tt=cputime;
load ui.dat; %load input
u=ui(:,1);
mx=5;
y(1:mx+1)=0;
for t=mx+1:500
for i=1:500
rand('twister',sum(100*clock));
e(i)=(rand-0.5)*0.002;
end
y(t)=0.797*y(t-1)-0.237*y(t-3)+0.441*u(t-1)+0.105*y(t-4)*u(t-4)...
+0.333*u(t-3)*u(t-5)+e(t);
% y(t)=1.5*y(t-1)-0.7*y(t-2)+u(t-5)+0.5*u(t-6);%+e(t);
% y(t)=0.5*y(t-1)-0.45*y(t-3)+1.2*u(t-3)+0.1*u(t-4);%+e(t); %Simulation 3
% y(t)=0.3*u(t-1)*u(t-2)-0.8*y(t-2)*u(t-1)*u(t-1)-0.1*u(t-1)*u(t-1)*u(t-2);%+e(t);
end
y=y';%load output
nodata=500;
noy=5; nou=5;
nl=2; start=max(nou,noy)+1;
m=nou+noy;
mphi=[]; %form regressor for ARMAX model
for i=start:nodata,
phi=[];
if nl==1
for j=1:noy,
phi=[phi y(i-j)];
end
for k=0:nou,
phi=[phi u(i-k)];
end;
end
if nl==2
for j=1:noy,
phi=[phi y(i-j)];
end
for k=0:nou,
phi=[phi u(i-k)];
end;
for j=1:noy,
for l=j:noy,
phi=[phi y(i-j)*y(i-l)];
end;
for k=1:nou,
phi=[phi y(i-j)*u(i-k)];
end;
end;
for k=1:nou,
for n=k:nou;
phi=[phi u(i-k)*u(i-n)];
end;
end;
end;
if nl==3
for j=1:noy,
phi=[phi y(i-j)];
end
for k=0:nou,
phi=[phi u(i-k)];
end;
for j=1:noy,
for l=j:noy,
phi=[phi y(i-j)*y(i-l)];
end;
for k=1:nou,
phi=[phi y(i-j)*u(i-k)];
end;
end;
for k=1:nou,
for n=k:nou;
phi=[phi u(i-k)*u(i-n)];
end;
end;
for j=1:noy,
for l=j:noy,
for p=l:noy
phi=[phi y(i-j)*y(i-l)*y(i-p)];
end;
for k=1:nou
phi=[phi y(i-j)*y(i-l)*u(i-k)];
end;
end;
for k=1:nou,
for n=k:nou,
phi=[phi y(i-j)*u(i-k)*u(i-n)];
end;
end;
end;
for k=1:nou,
for n=k:nou;
for q=n:nou
phi=[phi u(i-k)*u(i-n)*u(i-q)];
end
end;
end;
end;
mphi=[mphi;phi]; %end regressor
end;
[temp m]=size(mphi)
%% Initialization
%CONTROL PARAMETERS
clear
clc
n = 50; % Size of the swarm " no of birds "
bird_step = 50; % Maximum number of "birds steps"
dim = 10; % Dimension of the problem
c2 =1.2; % PSO parameter C1
c1 = 0.12; % PSO parameter C2
w =0.9; % pso momentum or inertia
fitness=0*ones(n,bird_step); %fitness population
%**********************************%
% Initialize the parameter
%**********************************%
R1 = rand(dim, n);
R2 = rand(dim, n);
current_fitness =0*ones(1,n);
%****************************************************%
% Initializing swarm and velocities and position
%****************************************************%
%Initialize random number of generator
current_position = 10*(rand(dim, n)-.5);
velocity = .3*randn(dim, n) ;
local_best_position = current_position ;
%**********************************%
% Evaluate initial population
%**********************************%
for i = 1:n
current_fitness(i) = rosenbrock(current_position(:,i));
end
local_best_fitness = current_fitness ;
[global_best_fitness,g] = min(local_best_fitness) ;
for i=1:n
globl_best_position(:,i) = local_best_position(:,g) ;
end
%**********************************%
% Optimization
%**********************************%
%------VELOCITY UPDATE------%
velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(globl_best_position-current_position));
%------SWARMUPDATE------%
current_position = current_position + velocity ;
%------EVALUATE NEW SWARM------%
global_best_fitness = inf;
%% MAIN LOOP %%
iter = 0 ; % Iterationscounter
while ( iter < bird_step )
iter = iter + 1;
for i = 1:n,
current_fitness(i) = rosenbrock(current_position(:,i)) ;
end
for i = 1 : n
if current_fitness(i) < local_best_fitness(i)
local_best_fitness(i) = current_fitness(i);
local_best_position(:,i) = current_position(:,i) ;
end
end
[current_global_best_fitness,g] = min(local_best_fitness);
if current_global_best_fitness < global_best_fitness
global_best_fitness = current_global_best_fitness;
for i=1:n
globl_best_position(:,i) = local_best_position(:,g);
end
end
velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(globl_best_position-current_position));
current_position = current_position + velocity;
x=current_position(1,:);
y=current_position(2,:);
clf
plot(x, y , 'h')
axis([-5 5 -5 5]);
pause(.3)
end % end of while loop its mean the end of all step that the birds move it
[Jbest_min,I] = min(current_fitness) % minimum fitness
current_position(:,I) % best solution
%
Below is the data input under name ui.dat
-3.7287749e-001
1.5462074e-001
3.6396491e-001
-2.2540157e-001
3.4016394e-001
-4.2924476e-001
-1.2120792e-001
-2.3183057e-001
-3.4707882e-001
1.3099958e-001
-1.8362518e-001
4.5911185e-001
-1.3260108e-003
2.3860852e-001
-4.8724433e-001
1.0535319e-001
7.6451321e-002
3.0737882e-001
1.5496806e-001
3.7822840e-001
4.0237293e-001
-3.4776714e-001
-3.0741991e-001
2.9097565e-001
-4.3929490e-001
-1.1017281e-001
-2.0003415e-001
2.3417991e-001
-3.9579075e-001
2.9257549e-001
2.8272894e-001
3.2397693e-002
-2.4664766e-001
-4.2904541e-001
1.2580333e-001
-4.7531853e-001
-4.3795751e-001
-3.7038802e-001
-4.9385928e-002
1.7233584e-001
3.5611083e-001
-1.5545757e-003
-4.5121550e-001
-1.8616764e-001
1.4163117e-001
2.8638664e-001
-2.1084951e-001
-2.1319893e-003
3.1843462e-001
9.5128293e-002
3.6424609e-002
-1.6912686e-001
-8.8309774e-002
2.9400624e-001
-1.5679178e-001
-3.7393264e-002
-1.3217633e-001
1.7957062e-001
6.7779032e-002
1.5177685e-001
-8.8862517e-003
-1.0154177e-001
-2.2515791e-002
-4.3341238e-001
-8.8971513e-002
4.6909196e-001
2.8072103e-001
2.2901789e-001
2.6565062e-001
2.5658080e-001
3.4327032e-001
2.7015685e-001
4.7866093e-001
-3.8864056e-001
-1.0392306e-001
-7.9393492e-003
-2.4190538e-001
-4.6303378e-001
4.7438132e-001
2.2642598e-001
-3.5203517e-001
-3.5210544e-001
2.0482410e-001
-1.1900648e-001
-4.2358856e-001
-8.9159657e-002
-3.5700945e-001
2.9892012e-001
4.3024816e-001
-4.9528052e-001
1.5004058e-001
1.7853246e-001
-2.4637716e-001
3.4318091e-001
-2.0604123e-001
-4.7313962e-001
-4.0669284e-001
2.9788497e-001
2.1140263e-001
2.8340678e-001
1.2392264e-001
3.2541546e-001
-4.6497741e-001
-9.4525028e-002
-2.5033124e-001
-1.9103158e-002
3.8083913e-001
-2.1931518e-001
9.9142509e-002
-4.7377386e-001
-3.4480133e-001
3.3391087e-001
-3.0510663e-001
3.2978804e-001
-1.6192514e-001
1.7111655e-001
-4.4763368e-001
2.3430907e-001
-5.1822218e-004
4.4329994e-001
-2.1022758e-001
-1.2343957e-001
-3.8624464e-001
4.6485071e-001
-6.7488433e-002
-4.1543680e-001
2.1667641e-001
6.7873795e-003
-1.7191226e-001
2.5352414e-001
3.3600582e-001
-2.4628365e-001
3.4425017e-002
-6.4834767e-002
-3.4229585e-001
1.0048120e-001
4.3745076e-001
-3.9224142e-001
3.9998143e-001
5.0464991e-002
-7.2643313e-002
-3.4761917e-001
-2.5245400e-001
-5.2631206e-002
3.2783438e-002
-1.4534887e-001
2.7311481e-001
3.8168114e-001
2.3409209e-001
-9.3567794e-002
1.0417875e-001
1.4110505e-001
-3.7253316e-001
-3.8080031e-003
-1.8953064e-001
7.8574464e-002
4.4360894e-001
-7.3060391e-002
-4.6687156e-001
4.2943388e-001
4.2498225e-001
-1.4169308e-001
-2.4001088e-001
2.8685965e-001
1.1582243e-002
6.2527153e-002
1.8479376e-001
-4.0760340e-001
3.7257897e-001
4.4293788e-001
-4.0340576e-001
3.4588251e-001
4.0939561e-001
-4.8865886e-001
2.3682704e-002
1.5034527e-001
-1.1485465e-001
1.4930205e-001
2.6285348e-001
7.5685981e-002
1.3191915e-001
-2.2179600e-001
3.3983642e-001
-7.3165044e-002
1.3162238e-001
3.3346698e-001
-2.2981454e-001
-9.9199264e-002
5.4254592e-002
-5.6134643e-002
-4.0961588e-001
2.4438131e-001
-4.6738455e-001
-7.0257457e-002
-4.6273772e-001
4.7579775e-001
2.2339908e-002
4.0962998e-001
-1.1675208e-001
3.8445224e-001
-2.4498239e-001
4.0904984e-001
3.9456044e-001
-1.0148284e-001
1.2502000e-001
6.7597286e-002
3.9451189e-001
-2.8583420e-001
-4.9614084e-001
3.8058146e-001
-2.6487836e-001
-2.5513538e-001
1.4091704e-001
-1.9547552e-001
3.2562081e-001
3.8368871e-001
4.4537296e-001
-1.0920466e-001
3.0132006e-001
-3.4288671e-001
1.2516540e-001
1.9898535e-001
-4.1413107e-001
3.1179573e-002
3.8856652e-001
-2.3633399e-001
-2.6521447e-001
3.3965723e-001
-4.4599703e-003
-3.4763363e-001
-2.6923184e-001
1.5795396e-001
6.2948381e-002
-2.0817126e-001
1.2230498e-001
2.1590541e-001
-2.1926897e-001
-8.7726841e-002
-1.3779429e-001
2.8139202e-001
-3.6451285e-001
4.0206951e-001
-2.1036391e-001
-4.4207389e-004
2.8359039e-001
1.7706206e-001
-3.5018827e-001
1.9661948e-001
-3.7098696e-001
4.4594629e-001
3.8641185e-001
1.5004613e-002
1.7940875e-001
4.7679336e-001
-3.7453964e-001
2.5224305e-001
3.2705191e-001
2.8143033e-001
-3.0911733e-001
-7.1358659e-002
-4.8554325e-001
-1.7471451e-001
-3.6529797e-001
-4.9482538e-002
7.2274930e-002
2.9202328e-001
-8.0263390e-002
3.2536544e-002
4.2570439e-001
3.9908184e-001
4.4833398e-002
4.0112395e-001
-4.4817329e-001
3.0860252e-001
-1.6509542e-001
-2.7131688e-001
3.2240181e-001
-1.5176542e-001
-3.3452940e-001
-4.7186643e-001
4.5536977e-001
1.8029071e-001
3.6056211e-001
4.3909352e-001
1.8019396e-001
4.1742375e-001
-2.4330833e-001
3.8561813e-001
4.2004279e-001
-1.9993656e-001
-4.2660905e-001
2.6739873e-001
-4.1504778e-001
2.2876408e-001
-5.2110356e-002
1.5124866e-001
-3.3049813e-001
3.1449241e-002
1.3380086e-001
-4.8590400e-001
-2.9628665e-002
3.8632597e-001
-3.8597156e-001
-5.7458848e-002
1.5954789e-001
-2.0522703e-001
4.5036840e-001
1.9428625e-001
-2.9319347e-001
5.4762326e-002
3.7927848e-001
5.7857691e-002
2.5233406e-001
3.9490079e-001
3.4184417e-001
-3.6914332e-001
-3.1084607e-001
-3.4636017e-001
-4.7109793e-001
-4.9091510e-001
9.6451742e-002
1.0904932e-001
4.1892343e-001
2.3357444e-001
-1.9885305e-001
-4.4255210e-003
-2.4183698e-001
2.3285436e-001
-3.8323925e-001
2.4604156e-001
3.0978911e-001
2.4523368e-001
-1.6285669e-001
8.4324980e-002
-3.1047895e-002
-4.1273497e-001
3.2871743e-001
1.8594471e-001
-2.3267494e-001
4.6948356e-001
-3.1622580e-001
-2.0005908e-001
-8.8815109e-002
-2.6351081e-001
-3.0494468e-001
2.0537928e-001
-3.1945175e-001
2.2332918e-002
-2.0382818e-001
-3.7218112e-002
4.2523242e-001
-2.8411086e-001
-4.9899112e-001
4.0660624e-001
1.8004195e-001
1.4951805e-002
2.2068881e-002
-3.9708142e-001
4.9688049e-001
-1.4103118e-001
1.2524260e-001
-1.0663348e-001
-4.9233755e-001
4.5284934e-002
9.1103744e-003
-2.5322108e-001
-4.5460473e-001
3.4172974e-001
-4.5175008e-001
-1.8367950e-001
2.8341943e-001
4.7240006e-001
8.6464081e-002
2.7804385e-001
2.2769744e-001
1.5098680e-001
1.6461502e-001
4.3877980e-001
3.5081266e-002
-1.0155988e-001
1.7045766e-001
-5.9465528e-002
-3.6712606e-001
-6.0796252e-002
4.7642742e-002
-1.0486362e-001
-1.0172801e-001
2.5134865e-001
2.2350346e-002
-9.5669663e-003
-4.1132107e-001
-2.4914837e-001
-5.2440977e-002
1.3796090e-001
2.0944529e-001
4.9261981e-001
4.3219487e-001
-4.0777050e-001
4.5353997e-001
-3.3720434e-001
4.7045509e-001
9.7006665e-002
-2.5977288e-001
-4.2970484e-001
-1.9995878e-001
3.1354482e-001
-4.2329013e-001
-1.4552690e-001
-3.6798913e-001
-3.4182037e-001
-4.3785295e-001
2.0184348e-001
-4.1351830e-001
1.1678724e-001
-3.2622842e-001
1.5140106e-001
-1.3035940e-003
-2.1548931e-001
3.3055974e-001
3.1835988e-001
4.3816940e-001
-4.9967390e-001
1.4038922e-001
-4.9264431e-001
-3.9357889e-001
-3.9320559e-001
-1.3289096e-001
-2.6039230e-001
-1.5385980e-001
-2.5038016e-001
-1.1293569e-001
-7.8961630e-002
1.4007729e-001
2.8755297e-001
-2.3000558e-001
3.4398236e-001
2.4046838e-001
3.2610195e-001
-3.1780775e-001
-4.3456382e-001
1.1035026e-001
2.0155281e-001
-3.8838229e-001
-4.0417627e-001
9.7833772e-002
3.1223261e-001
3.1457817e-001
-4.1056303e-001
2.3127871e-001
4.0385658e-001
-4.7768213e-002
-4.2931232e-001
-2.5872184e-001
2.3186506e-001
-4.5950838e-001
-7.5477237e-002
4.0215087e-002
4.5382784e-001
-2.9109426e-001
-3.8366790e-001
1.4622007e-001
-3.9158891e-001
4.8349692e-001
-2.5165596e-001
1.0635572e-001
3.1669527e-001
3.3005421e-001
-1.0960668e-002
2.6072826e-001
4.1510794e-001
4.0097496e-001
-2.8576249e-001
4.7061106e-002
2.8470948e-001
-3.0555891e-001
2.4688915e-001
-2.4441776e-002
8.3258581e-002
-2.3945192e-001
-4.1517754e-001
-2.0186715e-001
4.1712936e-001
-2.9481807e-002
-2.3053245e-001
2.6297048e-001
2.7217211e-001
-4.7870038e-001
3.7998619e-001
2.9817018e-001
-1.7583476e-001
1.6904442e-001
-2.0370614e-001
4.2995204e-001
-2.1803957e-001
-3.3112003e-001
2.4516660e-001
-2.2865648e-002
1.5344454e-001
4.6657593e-001
-1.8697251e-001

Answers (1)

Sam Chak
Sam Chak on 23 Jun 2022
Your code is quite long. By the way, you are using Rosenbrock function as your fitness function, not your NARX model.
for i = 1:n
current_fitness(i) = rosenbrock(current_position(:,i));
end
  5 Comments
Paul
Paul on 23 Jun 2022
I'm afraid I can't help with this question.
I did notice this line
current_position = current_position + velocity;
which reads a bit odd. Perhaps "velocity" in this context means "change in position"?
Sam Chak
Sam Chak on 23 Jun 2022
No worries, @Paul. Thanks for checking out the question.
Regarding this:
current_position = current_position + velocity;
I have checked the particle swarm algorithm in MATLAB doc:
Another source also says the same thing:
I guess the update theory is similar to the Calculus or Finite Difference:
dx/dt = v
dx = v*dt
dx = x2 - x1
x2 = x1 + v*dt

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!