clc
clear
close all
%% ===============================
% READ HEAT EXCHANGER DATASET
%% ===============================
data = readmatrix(‘exchanger.xlsx’);
% Column 1 = Input
g = data(:,1);
% Column 2 = Output
x = data(:,2);
% Number of estimation samples
M = 3863;
% Select first M samples
g = g(1:M);
x = x(1:M);
% Convert to column vectors
g = g(:);
x = x(:);
%% ===============================
% PSO PARAMETERS
%% ===============================
dim = 8;
nParticles = 100;
maxIter = 800;
lb = [-2.0 0.3 0.0 -0.4 -40 30 -30 55];
ub = [-1.0 0.9 0.4 -0.1 -25 50 -20 75];
%% ===============================
% INITIALIZATION
%% ===============================
particle = struct;
for i = 1:nParticles
particle(i).position = lb + rand(1,dim).*(ub-lb);
particle(i).velocity = zeros(1,dim);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
costs = [particle.cost];
[~,idx] = min(costs);
global_best = particle(idx).best;
%% ===============================
% PSO PARAMETERS
%% ===============================
w = 0.6;
c1 = 2.5;
c2 = 2.5;
BestCost = zeros(maxIter,1);
theta_history = zeros(maxIter,dim);
mse_history = zeros(maxIter,1);
%% ===============================
% PSO MAIN LOOP
%% ===============================
for iter = 1:maxIter
for i = 1:nParticles
particle(i).velocity = …
w*particle(i).velocity …
+ c1*rand*(particle(i).best.position-particle(i).position) …
+ c2*rand*(global_best.position-particle(i).position);
particle(i).position = …
particle(i).position + particle(i).velocity;
particle(i).position = …
max(particle(i).position,lb);
particle(i).position = …
min(particle(i).position,ub);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
if particle(i).cost < particle(i).best.cost
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
if particle(i).best.cost < global_best.cost
global_best = particle(i).best;
end
end
theta_history(iter,:) = global_best.position;
mse_history(iter) = global_best.cost;
BestCost(iter) = global_best.cost;
fprintf(‘Iteration %d Best Error = %en’, …
iter,global_best.cost);
end
%% ===============================
% Estimated Parameters
%% ===============================
xext = global_best.position;
disp(‘Estimated Parameters’)
disp(xext)
true_theta = [-1.6080 0.6385 0.2070 -0.1764 -31.55 41.73 -24.20 68.03];
disp(‘True Parameters’)
disp(true_theta)
%% ===============================
% Parameter Extraction
%% ===============================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ===============================
% Nonlinear Block
%% ===============================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ===============================
% Estimated Output
%% ===============================
xhat = zeros(M,1);
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1)-0.3*xhat(t-2)-0.105*xhat(t-3)-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2)-0.3*xhat(t-3)-0.105*xhat(t-4)-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1)-0.3*khat(t-2)-0.105*khat(t-3)-0.0595*khat(t-4)) …
+ delta2*(khat(t-2)-0.3*khat(t-3)-0.105*khat(t-4)-0.0595*khat(t-5));
end
%% ===============================
% Performance
%% ===============================
disp(‘Best Cost from PSO’)
disp(global_best.cost)
err = x – xhat;
valid = isfinite(err) & (abs(err) < 1e6);
if any(valid)
mse_value = mean(err(valid).^2);
else
mse_value = Inf;
end
fprintf(‘Final MSE = %en’, mse_value);
%% ===============================
% Actual vs Estimated Output
%% ===============================
figure
plot(x,’b’,’LineWidth’,1.5)
hold on
plot(xhat,’r–‘,’LineWidth’,1.5)
xlabel(‘Samples’)
ylabel(‘Output’)
title(‘Actual vs Estimated Output’)
legend(‘Actual Output’,’Estimated Output’)
grid on
function j = fitness_function(xext,x,g,M)
%% ==========================================
% Column vectors
%% ==========================================
x = x(:);
g = g(:);
%% ==========================================
% Parameters
%% ==========================================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ==========================================
% Nonlinear Block
%% ==========================================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ==========================================
% Free Run Simulation
%% ==========================================
xhat = zeros(M,1);
% Initial conditions
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1) …
-0.3*xhat(t-2) …
-0.105*xhat(t-3) …
-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2) …
-0.3*xhat(t-3) …
-0.105*xhat(t-4) …
-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1) …
-0.3*khat(t-2) …
-0.105*khat(t-3) …
-0.0595*khat(t-4)) …
+ delta2*(khat(t-2) …
-0.3*khat(t-3) …
-0.105*khat(t-4) …
-0.0595*khat(t-5));
end
%% ==========================================
% Mean Squared Error
%% ==========================================
err = x – xhat;
% Keep only reasonable errors
valid = isfinite(err) & (abs(err) < 1e6);
if sum(valid) < 0.80*M
j = 1e20;
else
err = err(valid);
j = mean(err.^2);
end
if isnan(j) || isinf(j)
j = 1e20;
end
endclc
clear
close all
%% ===============================
% READ HEAT EXCHANGER DATASET
%% ===============================
data = readmatrix(‘exchanger.xlsx’);
% Column 1 = Input
g = data(:,1);
% Column 2 = Output
x = data(:,2);
% Number of estimation samples
M = 3863;
% Select first M samples
g = g(1:M);
x = x(1:M);
% Convert to column vectors
g = g(:);
x = x(:);
%% ===============================
% PSO PARAMETERS
%% ===============================
dim = 8;
nParticles = 100;
maxIter = 800;
lb = [-2.0 0.3 0.0 -0.4 -40 30 -30 55];
ub = [-1.0 0.9 0.4 -0.1 -25 50 -20 75];
%% ===============================
% INITIALIZATION
%% ===============================
particle = struct;
for i = 1:nParticles
particle(i).position = lb + rand(1,dim).*(ub-lb);
particle(i).velocity = zeros(1,dim);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
costs = [particle.cost];
[~,idx] = min(costs);
global_best = particle(idx).best;
%% ===============================
% PSO PARAMETERS
%% ===============================
w = 0.6;
c1 = 2.5;
c2 = 2.5;
BestCost = zeros(maxIter,1);
theta_history = zeros(maxIter,dim);
mse_history = zeros(maxIter,1);
%% ===============================
% PSO MAIN LOOP
%% ===============================
for iter = 1:maxIter
for i = 1:nParticles
particle(i).velocity = …
w*particle(i).velocity …
+ c1*rand*(particle(i).best.position-particle(i).position) …
+ c2*rand*(global_best.position-particle(i).position);
particle(i).position = …
particle(i).position + particle(i).velocity;
particle(i).position = …
max(particle(i).position,lb);
particle(i).position = …
min(particle(i).position,ub);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
if particle(i).cost < particle(i).best.cost
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
if particle(i).best.cost < global_best.cost
global_best = particle(i).best;
end
end
theta_history(iter,:) = global_best.position;
mse_history(iter) = global_best.cost;
BestCost(iter) = global_best.cost;
fprintf(‘Iteration %d Best Error = %en’, …
iter,global_best.cost);
end
%% ===============================
% Estimated Parameters
%% ===============================
xext = global_best.position;
disp(‘Estimated Parameters’)
disp(xext)
true_theta = [-1.6080 0.6385 0.2070 -0.1764 -31.55 41.73 -24.20 68.03];
disp(‘True Parameters’)
disp(true_theta)
%% ===============================
% Parameter Extraction
%% ===============================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ===============================
% Nonlinear Block
%% ===============================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ===============================
% Estimated Output
%% ===============================
xhat = zeros(M,1);
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1)-0.3*xhat(t-2)-0.105*xhat(t-3)-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2)-0.3*xhat(t-3)-0.105*xhat(t-4)-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1)-0.3*khat(t-2)-0.105*khat(t-3)-0.0595*khat(t-4)) …
+ delta2*(khat(t-2)-0.3*khat(t-3)-0.105*khat(t-4)-0.0595*khat(t-5));
end
%% ===============================
% Performance
%% ===============================
disp(‘Best Cost from PSO’)
disp(global_best.cost)
err = x – xhat;
valid = isfinite(err) & (abs(err) < 1e6);
if any(valid)
mse_value = mean(err(valid).^2);
else
mse_value = Inf;
end
fprintf(‘Final MSE = %en’, mse_value);
%% ===============================
% Actual vs Estimated Output
%% ===============================
figure
plot(x,’b’,’LineWidth’,1.5)
hold on
plot(xhat,’r–‘,’LineWidth’,1.5)
xlabel(‘Samples’)
ylabel(‘Output’)
title(‘Actual vs Estimated Output’)
legend(‘Actual Output’,’Estimated Output’)
grid on
function j = fitness_function(xext,x,g,M)
%% ==========================================
% Column vectors
%% ==========================================
x = x(:);
g = g(:);
%% ==========================================
% Parameters
%% ==========================================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ==========================================
% Nonlinear Block
%% ==========================================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ==========================================
% Free Run Simulation
%% ==========================================
xhat = zeros(M,1);
% Initial conditions
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1) …
-0.3*xhat(t-2) …
-0.105*xhat(t-3) …
-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2) …
-0.3*xhat(t-3) …
-0.105*xhat(t-4) …
-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1) …
-0.3*khat(t-2) …
-0.105*khat(t-3) …
-0.0595*khat(t-4)) …
+ delta2*(khat(t-2) …
-0.3*khat(t-3) …
-0.105*khat(t-4) …
-0.0595*khat(t-5));
end
%% ==========================================
% Mean Squared Error
%% ==========================================
err = x – xhat;
% Keep only reasonable errors
valid = isfinite(err) & (abs(err) < 1e6);
if sum(valid) < 0.80*M
j = 1e20;
else
err = err(valid);
j = mean(err.^2);
end
if isnan(j) || isinf(j)
j = 1e20;
end
end clc
clear
close all
%% ===============================
% READ HEAT EXCHANGER DATASET
%% ===============================
data = readmatrix(‘exchanger.xlsx’);
% Column 1 = Input
g = data(:,1);
% Column 2 = Output
x = data(:,2);
% Number of estimation samples
M = 3863;
% Select first M samples
g = g(1:M);
x = x(1:M);
% Convert to column vectors
g = g(:);
x = x(:);
%% ===============================
% PSO PARAMETERS
%% ===============================
dim = 8;
nParticles = 100;
maxIter = 800;
lb = [-2.0 0.3 0.0 -0.4 -40 30 -30 55];
ub = [-1.0 0.9 0.4 -0.1 -25 50 -20 75];
%% ===============================
% INITIALIZATION
%% ===============================
particle = struct;
for i = 1:nParticles
particle(i).position = lb + rand(1,dim).*(ub-lb);
particle(i).velocity = zeros(1,dim);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
costs = [particle.cost];
[~,idx] = min(costs);
global_best = particle(idx).best;
%% ===============================
% PSO PARAMETERS
%% ===============================
w = 0.6;
c1 = 2.5;
c2 = 2.5;
BestCost = zeros(maxIter,1);
theta_history = zeros(maxIter,dim);
mse_history = zeros(maxIter,1);
%% ===============================
% PSO MAIN LOOP
%% ===============================
for iter = 1:maxIter
for i = 1:nParticles
particle(i).velocity = …
w*particle(i).velocity …
+ c1*rand*(particle(i).best.position-particle(i).position) …
+ c2*rand*(global_best.position-particle(i).position);
particle(i).position = …
particle(i).position + particle(i).velocity;
particle(i).position = …
max(particle(i).position,lb);
particle(i).position = …
min(particle(i).position,ub);
particle(i).cost = fitness_function( …
particle(i).position,x,g,M);
if particle(i).cost < particle(i).best.cost
particle(i).best.position = particle(i).position;
particle(i).best.cost = particle(i).cost;
end
if particle(i).best.cost < global_best.cost
global_best = particle(i).best;
end
end
theta_history(iter,:) = global_best.position;
mse_history(iter) = global_best.cost;
BestCost(iter) = global_best.cost;
fprintf(‘Iteration %d Best Error = %en’, …
iter,global_best.cost);
end
%% ===============================
% Estimated Parameters
%% ===============================
xext = global_best.position;
disp(‘Estimated Parameters’)
disp(xext)
true_theta = [-1.6080 0.6385 0.2070 -0.1764 -31.55 41.73 -24.20 68.03];
disp(‘True Parameters’)
disp(true_theta)
%% ===============================
% Parameter Extraction
%% ===============================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ===============================
% Nonlinear Block
%% ===============================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ===============================
% Estimated Output
%% ===============================
xhat = zeros(M,1);
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1)-0.3*xhat(t-2)-0.105*xhat(t-3)-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2)-0.3*xhat(t-3)-0.105*xhat(t-4)-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1)-0.3*khat(t-2)-0.105*khat(t-3)-0.0595*khat(t-4)) …
+ delta2*(khat(t-2)-0.3*khat(t-3)-0.105*khat(t-4)-0.0595*khat(t-5));
end
%% ===============================
% Performance
%% ===============================
disp(‘Best Cost from PSO’)
disp(global_best.cost)
err = x – xhat;
valid = isfinite(err) & (abs(err) < 1e6);
if any(valid)
mse_value = mean(err(valid).^2);
else
mse_value = Inf;
end
fprintf(‘Final MSE = %en’, mse_value);
%% ===============================
% Actual vs Estimated Output
%% ===============================
figure
plot(x,’b’,’LineWidth’,1.5)
hold on
plot(xhat,’r–‘,’LineWidth’,1.5)
xlabel(‘Samples’)
ylabel(‘Output’)
title(‘Actual vs Estimated Output’)
legend(‘Actual Output’,’Estimated Output’)
grid on
function j = fitness_function(xext,x,g,M)
%% ==========================================
% Column vectors
%% ==========================================
x = x(:);
g = g(:);
%% ==========================================
% Parameters
%% ==========================================
sigma1 = xext(1);
sigma2 = xext(2);
delta1 = xext(3);
delta2 = xext(4);
l1 = xext(5);
l2 = xext(6);
l3 = xext(7);
l4 = xext(8);
%% ==========================================
% Nonlinear Block
%% ==========================================
khat = l1*g + l2*g.^2 + l3*g.^3 + l4*g.^4;
%% ==========================================
% Free Run Simulation
%% ==========================================
xhat = zeros(M,1);
% Initial conditions
xhat(1:5) = x(1:5);
for t = 6:M
xhat(t) = …
sigma1*(xhat(t-1) …
-0.3*xhat(t-2) …
-0.105*xhat(t-3) …
-0.0595*xhat(t-4)) …
+ sigma2*(xhat(t-2) …
-0.3*xhat(t-3) …
-0.105*xhat(t-4) …
-0.0595*xhat(t-5)) …
+ delta1*(khat(t-1) …
-0.3*khat(t-2) …
-0.105*khat(t-3) …
-0.0595*khat(t-4)) …
+ delta2*(khat(t-2) …
-0.3*khat(t-3) …
-0.105*khat(t-4) …
-0.0595*khat(t-5));
end
%% ==========================================
% Mean Squared Error
%% ==========================================
err = x – xhat;
% Keep only reasonable errors
valid = isfinite(err) & (abs(err) < 1e6);
if sum(valid) < 0.80*M
j = 1e20;
else
err = err(valid);
j = mean(err.^2);
end
if isnan(j) || isinf(j)
j = 1e20;
end
end is my code correct can anyone help me MATLAB Answers — New Questions
