Error using autoclave function
Show older comments
Hi, I am trying to run a simple particle filter function, but the graph that I get is not smooth like I expected. I tred to change the process equation to autoclave as I saw a good particle filter function using process equation autoclave. But there is an error when I run it. The error is
Arrays have incompatible sizes for this operation.
Error in autoclave (line 46)
xkm = [h Vf']+uk';
Error in testing8 (line 44)
x(:,k) = sys(k, x(:,k-1), u(:,k)); % simulate state
Any recommendation to solve this problem is really appreciated.
%% clear memory, screen, and close all figures
clear, clc, close all;
%% Process equation x[k] = sys(k, x[k-1], u[k]);% untuk state vector
nx = 50; % number of states
sys= @autoclave; % (returns column vector)
%% Observation equation y[k] = obs(k, x[k], v[k]);
ny = 1; % number of observations
obs = @(k, xk,vk) xk(1) + vk; % (returns column vector)
%% PDF of process noise and noise generator function
nu = 50; % size of the vector of process noise
sigma_u = sqrt(10);
p_sys_noise = @(u) normpdf(u, 0, sigma_u);
gen_sys_noise = @(u) normrnd(0, sigma_u); % sample from p_sys_noise (returns column vector)
%% PDF of observation noise and noise generator function
nv = 1; % size of the vector of observation noise
sigma_v = sqrt(10);
p_obs_noise = @(v) normpdf(v, 0, sigma_v);
gen_obs_noise = @(v) normrnd(0, sigma_v); % sample from p_obs_noise (returns column vector)
%% Initial PDF
% p_x0 = @(x) normpdf(x, 0,sqrt(10)); % initial pdf
gen_x0 = @(x) [89.4;93.75;90;94.38;98.75;94.38;96.25;94.38;96.25;100;95;96.25;95;95;98.75;98.125;97.5;98.75;100;100;91.875;93.75;93.125;96.25;98.125;98.125;99.375;100;100;100;93.75;92.5;96.25;95.625;99.375;97.5;100;99.375;99.375;100;95.625;95.625;97.5;97.5;98.75;99.375;99.375;91.25;98.125;100]+ones(50,1)*normrnd(0, sqrt(10)); % sample from p_x0 (returns column vector)
%% Transition prior PDF p(x[k] | x[k-1])
% (under the suposition of additive process noise)
% p_xk_given_xkm1 = @(k, xk, xkm1) p_sys_noise(xk - sys(k, xkm1, 0));
%% Observation likelihood PDF p(y[k] | x[k])
% (under the suposition of additive process noise)
p_yk_given_xk = @(k, yk, xk) p_obs_noise(yk - obs(k, xk, 0));
%% %% Number of time steps
T = 100;
%% Separate memory space
x = zeros(nx,T); y = zeros(ny,T);
u = zeros(nu,T); v = zeros(nv,T);
%% Simulate system
xh0 = [89.4;93.75;90;94.38;98.75;94.38;96.25;94.38;96.25;100;95;96.25;95;95;98.75;98.125;97.5;98.75;100;100;91.875;93.75;93.125;96.25;98.125;98.125;99.375;100;100;100;93.75;92.5;96.25;95.625;99.375;97.5;100;99.375;99.375;100;95.625;95.625;97.5;97.5;98.75;99.375;99.375;91.25;98.125;100]; % initial state
u(:,1) = 0; % initial process noise
v(:,1) = gen_obs_noise(sigma_v); % initial observation noise
x(:,1) = xh0;
y(:,1) = obs(1, xh0, v(:,1));
for k = 2:T
% here we are basically sampling from p_xk_given_xkm1 and from p_yk_given_xk
u(:,k) = 0.1*gen_sys_noise(); % simulate process noise
v(:,k) = gen_obs_noise(); % simulate observation noise
x(:,k) = sys(k, x(:,k-1), u(:,k)); % simulate state
y(:,k) = obs(k, x(:,k), v(:,k)); % simulate observation
end
fprintf('Finish simulate system \n')
%% Separate memory
xh = zeros(nx, T); xh(:,1) = xh0;
yh = zeros(ny, T); yh(:,1) = obs(1, xh0, 0);
pf.k = 1; % initial iteration number
pf.Ns = 100; % number of particles
pf.w = zeros(pf.Ns, T); % weights
pf.particles = zeros(nx, pf.Ns, T); % particles
pf.gen_x0 = gen_x0; % function for sampling from initial pdf p_x0
pf.p_yk_given_xk = p_yk_given_xk; % function of the observation likelihood PDF p(y[k] | x[k])
pf.gen_sys_noise = gen_sys_noise; % function for generating system noise
%pf.p_x0 = p_x0; % initial prior PDF p(x[0])
%pf.p_xk_given_ xkm1 = p_xk_given_xkm1; % transition prior PDF p(x[k] | x[k-1])
%% Estimate state
for k = 2:T
fprintf('Iteration = %d/%d\n',k,T);
% state estimation
pf.k = k;
%[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
% filtered observation
yh(:,k) = obs(k, xh(:,k), 0);
end
%% Make plots of the evolution of the density
figure
hold on;
xi = 1:T;
yi = -25:0.25:25;
[xx,yy] = meshgrid(xi,yi);
den = zeros(size(xx));
xhmode = zeros(size(xh));
for i = xi
% for each time step perform a kernel density estimation
den(:,i) = ksdensity(pf.particles(1,:,i), yi,'kernel','epanechnikov');
[~, idx] = max(den(:,i));
% estimate the mode of the density
xhmode(i) = yi(idx);
plot3(repmat(xi(i),length(yi),1), yi', den(:,i));
end
view(3);
box on;
title('Evolution of the state density','FontSize',14)
figure
mesh(xx,yy,den);
title('Evolution of the state density','FontSize',14)
%% plot of the state vs estimated state by the particle filter vs particle paths
figure
hold on;
%h1 = plot(1:T,squeeze(pf.particles),'y');
h1 = plot(1:T,squeeze(pf.particles(2,:,:)),'y');
h2 = plot(1:T,x(1,:),'b','LineWidth',1);
h3 = plot(1:T,xh(1,:),'r','LineWidth',1);
h4 = plot(1:T,y(1,:),'g.','LineWidth',1);
legend([h2 h3 h4 h1(1)],'state','mean of estimated state','mode of estimated state','particle paths');
title('State vs estimated state by the particle filter vs particle paths','FontSize',14);
%% plot of the observation vs filtered observation by the particle filter
figure
plot(1:T,y,'b', 1:T,yh,'r');
legend('observation','filtered observation');
title('Observation vs filtered observation by the particle filter','FontSize',14);
return;
4 Comments
Image Analyst
on 22 Nov 2022
You forgot to attach your autoclave function.
シティヌルシュハダ モハマド ナシル
on 23 Nov 2022
Image Analyst
on 23 Nov 2022
I don't have it and it doesn't seem to be a built-in function as you can see by the red error text in your original post.
Of course I did not intend to attach it -- I have absolutely no idea what it is. I assume that you have the function, not us.
I do not see it in your code that you posted. So I don't know what it is, or where it is.
シティヌルシュハダ モハマド ナシル
on 24 Nov 2022
Edited: Walter Roberson
on 24 Nov 2022
Accepted Answer
Walter Roberson
on 24 Nov 2022
Moved: Walter Roberson
on 24 Nov 2022
%% clear memory, screen, and close all figures
clear, clc, close all;
%% Process equation x[k] = sys(k, x[k-1], u[k]);% untuk state vector
nx = 50; % number of states
sys= @autoclave; % (returns column vector)
%% Observation equation y[k] = obs(k, x[k], v[k]);
ny = 1; % number of observations
obs = @(k, xk,vk) xk(1) + vk; % (returns column vector)
%% PDF of process noise and noise generator function
nu = 50; % size of the vector of process noise
sigma_u = sqrt(10);
p_sys_noise = @(u) normpdf(u, 0, sigma_u);
gen_sys_noise = @(u) normrnd(0, sigma_u); % sample from p_sys_noise (returns column vector)
%% PDF of observation noise and noise generator function
nv = 1; % size of the vector of observation noise
sigma_v = sqrt(10);
p_obs_noise = @(v) normpdf(v, 0, sigma_v);
gen_obs_noise = @(v) normrnd(0, sigma_v); % sample from p_obs_noise (returns column vector)
%% Initial PDF
% p_x0 = @(x) normpdf(x, 0,sqrt(10)); % initial pdf
gen_x0 = @(x) [89.4;93.75;90;94.38;98.75;94.38;96.25;94.38;96.25;100;95;96.25;95;95;98.75;98.125;97.5;98.75;100;100;91.875;93.75;93.125;96.25;98.125;98.125;99.375;100;100;100;93.75;92.5;96.25;95.625;99.375;97.5;100;99.375;99.375;100;95.625;95.625;97.5;97.5;98.75;99.375;99.375;91.25;98.125;100]+ones(50,1)*normrnd(0, sqrt(10)); % sample from p_x0 (returns column vector)
%% Transition prior PDF p(x[k] | x[k-1])
% (under the suposition of additive process noise)
% p_xk_given_xkm1 = @(k, xk, xkm1) p_sys_noise(xk - sys(k, xkm1, 0));
%% Observation likelihood PDF p(y[k] | x[k])
% (under the suposition of additive process noise)
p_yk_given_xk = @(k, yk, xk) p_obs_noise(yk - obs(k, xk, 0));
%% %% Number of time steps
T = 100;
%% Separate memory space
x = zeros(nx,T); y = zeros(ny,T);
u = zeros(nu,T); v = zeros(nv,T);
%% Simulate system
xh0 = [89.4;93.75;90;94.38;98.75;94.38;96.25;94.38;96.25;100;95;96.25;95;95;98.75;98.125;97.5;98.75;100;100;91.875;93.75;93.125;96.25;98.125;98.125;99.375;100;100;100;93.75;92.5;96.25;95.625;99.375;97.5;100;99.375;99.375;100;95.625;95.625;97.5;97.5;98.75;99.375;99.375;91.25;98.125;100]; % initial state
u(:,1) = 0; % initial process noise
v(:,1) = gen_obs_noise(sigma_v); % initial observation noise
x(:,1) = xh0;
y(:,1) = obs(1, xh0, v(:,1));
for k = 2:T
% here we are basically sampling from p_xk_given_xkm1 and from p_yk_given_xk
u(:,k) = 0.1*gen_sys_noise(); % simulate process noise
v(:,k) = gen_obs_noise(); % simulate observation noise
x(:,k) = sys(k, x(:,k-1), u(:,k)); % simulate state
y(:,k) = obs(k, x(:,k), v(:,k)); % simulate observation
end
You can see that h is scalar, Vf is 5 x 1 so Vf' is 1 x 5, so [h Vf'] would be 1 x 6.
uk is 50 x 1, so uk' is 1 x 50.
You cannot add a 1 x 6 array and a 1 x 50 array.
The 1 x 5 for Vf is because Nz is 5.
fprintf('Finish simulate system \n')
%% Separate memory
xh = zeros(nx, T); xh(:,1) = xh0;
yh = zeros(ny, T); yh(:,1) = obs(1, xh0, 0);
pf.k = 1; % initial iteration number
pf.Ns = 100; % number of particles
pf.w = zeros(pf.Ns, T); % weights
pf.particles = zeros(nx, pf.Ns, T); % particles
pf.gen_x0 = gen_x0; % function for sampling from initial pdf p_x0
pf.p_yk_given_xk = p_yk_given_xk; % function of the observation likelihood PDF p(y[k] | x[k])
pf.gen_sys_noise = gen_sys_noise; % function for generating system noise
%pf.p_x0 = p_x0; % initial prior PDF p(x[0])
%pf.p_xk_given_ xkm1 = p_xk_given_xkm1; % transition prior PDF p(x[k] | x[k-1])
%% Estimate state
for k = 2:T
fprintf('Iteration = %d/%d\n',k,T);
% state estimation
pf.k = k;
%[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
[xh(:,k), pf] = particle_filter(sys, y(:,k), pf, 'multinomial_resampling');
% filtered observation
yh(:,k) = obs(k, xh(:,k), 0);
end
%% Make plots of the evolution of the density
figure
hold on;
xi = 1:T;
yi = -25:0.25:25;
[xx,yy] = meshgrid(xi,yi);
den = zeros(size(xx));
xhmode = zeros(size(xh));
for i = xi
% for each time step perform a kernel density estimation
den(:,i) = ksdensity(pf.particles(1,:,i), yi,'kernel','epanechnikov');
[~, idx] = max(den(:,i));
% estimate the mode of the density
xhmode(i) = yi(idx);
plot3(repmat(xi(i),length(yi),1), yi', den(:,i));
end
view(3);
box on;
title('Evolution of the state density','FontSize',14)
figure
mesh(xx,yy,den);
title('Evolution of the state density','FontSize',14)
%% plot of the state vs estimated state by the particle filter vs particle paths
figure
hold on;
%h1 = plot(1:T,squeeze(pf.particles),'y');
h1 = plot(1:T,squeeze(pf.particles(2,:,:)),'y');
h2 = plot(1:T,x(1,:),'b','LineWidth',1);
h3 = plot(1:T,xh(1,:),'r','LineWidth',1);
h4 = plot(1:T,y(1,:),'g.','LineWidth',1);
legend([h2 h3 h4 h1(1)],'state','mean of estimated state','mode of estimated state','particle paths');
title('State vs estimated state by the particle filter vs particle paths','FontSize',14);
%% plot of the observation vs filtered observation by the particle filter
figure
plot(1:T,y,'b', 1:T,yh,'r');
legend('observation','filtered observation');
title('Observation vs filtered observation by the particle filter','FontSize',14);
return;
function [xkm] = autoclave(k, xkm1, uk)
Vfo = 0.30; % initial fiber volume fraction
h0 = 0.02; % initial part thickness, [m]
mu = 0.40; % resin viscosity, [Pa.s]
C = 5.0e-10; % coefficient [m^2] in Kozeny-Carman permeability:
% Kzz = C (1-Vf)^3/Vf^2
Va = 0.70; % maximum attainable fiber volume fraction
As = 280; % coefficient [Pa] in Fiber in a Box model:
% Sigma_zz = As (sqrt(Vf/Vo)-1)/(sqrt(Va/Vf)-1)^4
Nz = 5; % # of nodes in z'
z = linspace(0,1,Nz); % z' array, [m]
dz = z(2)-z(1); % Delta z', [m]
dt=0.001;
h=xkm1(1);
Vf = xkm1(2:Nz+1); % fiber volume fraction
PT = 10.0e5; % autoclave pressure, [Pa]
e = (1-Vf)./Vf;
e_next = e;
e0 = e(1);
Kzz = C*(1-Vf).^3./Vf.^2; % permeability, [m^2]
Sigmazz = As*(sqrt(Vf./Vfo)-1)./(sqrt(Va./Vf)-1).^4;
VF = Vfo:0.0001:Va-0.0001;
SIGMAZZ = As*(sqrt(VF./Vfo)-1)./(sqrt(Va./VF)-1).^4;
e_next(1) = (1/3)*(-e(3)+4*e(2)); % B.C. at z = 0
c = VF(find(abs(SIGMAZZ-PT) == min(abs(SIGMAZZ-PT))));
e_next(end) = (1-c)/c; % B.C. at z = h
for i = 2:Nz-1
c1 = - ( (Kzz(i+1)/(1+e(i+1))) - (Kzz(i-1)/(1+e(i-1))))/(2*dz );
c2 = ( Sigmazz(i+1) - Sigmazz(i-1) )/(2*dz );
c3 = - Kzz(i)/(1+e(i));
c4 = ( Sigmazz(i+1) - 2*Sigmazz(i) + Sigmazz(i-1) )/( dz^2);
e_next(i) = e(i) + (dt*(1+e0^2)/mu) * (1/h^2) * (c1*c2+c3*c4);
end
Vf = 1./(1+e_next);
h = (h0*Vfo) * (1/mean(Vf));
whos
xkm = [h Vf']+uk';
end
1 Comment
シティヌルシュハダ モハマド ナシル
on 2 Dec 2022
More Answers (0)
Categories
Find more on Quadratic Programming and Cone Programming in Help Center and File Exchange
on 22 Nov 2022
on 2 Dec 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
