How to compute a model (Bouc-Wen) including a differential equation ?

Hello everyone,

I am having an hard time trying to build a Bouc-Wen model on MATLAB. Although it is easy to build on Simulink, I need to learn how to write it on a .m file in order to use Genetic Algorithms on it in the future.

The Bouc-Wen model is defined by the two equations below, one giving the force F in function of displacement and velocity x and dx, constant parameters c0, k, alpha, f0 and a variable z which is defined by a differential equation including velocity dx and constant parameters gamma, nu, n and A.

with :

In a first time, I am trying to compute the force F predicted by the model, for given sinusoidal displacement and velocity, and given constant parameters. To do this, I am trying in a first time to solve the differential equation then trying to use the solution to compute Fmodel

Here is my code using dsolve to solve the differential equation and computing Fmodel:

% parameters of sinusoidal displacement and velocity
Ampl = 0.013;
freqHz = 0.6;
t=0:1e-3:20;

% vectors describing virtual displacement and virtual velocity
X_virt = Ampl*sin(2*pi*freqHz*t);
V_virt = Ampl*2*pi*freqHz*cos(2*pi*freqHz*t); 

% vector of constant parameters
vp = [650 300 8e4 0 1.2e-6 1e-6 15 2];

% solving differential equation
dx = V_virt;
syms z(dx)
z(dx) = dsolve (diff(z) == -vp(5)*abs(dx)*z*(abs(z))^(vp(8)-1)...
            -vp(6)*dx*(abs(z))^(vp(8)) + vp(7)*dx)

 % computing FModel       
   Fmodel = vp(1)*V_virt + vp(2)*X_virt + vp(3)*z(dx) + vp(4);

% plotting
 figure('Name', 'Test BW');
 grid on;
 plot (V_virt, Fmodel)

And here is the error message I get :

Error using symengine (line 59)
Array sizes must match.

Error in sym/privBinaryOp (line 903)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  +  (line 7)
X = privBinaryOp(A, B, 'symobj::zip', '_plus');

Error in test_bw (line 20)
 Fmodel = vp(1)*V_virt + vp(2)*X_virt + vp(3)*z(dx) + vp(4);

The error may come from the fact that z(dx) is a sym variable so I can't use it like this. As I am a beginner on how to manipulate differential equations, any hints on how to proceed would be awesome.

Many thanks,

Pierre,

PS: I am using R2015a release.

[ain dolkifl]

Hello,

did you find how to compute a model (Bouc-Wen) including a differential equation in matlab beacause I am stuck at the same stage?

[ketan sengar]

If u have matlab code for bouc-wen model then please share with me. I am trying to build this code, but still not getting result.

Answer 1

[Wilson Rocha Lacerda Junior]

Take a look at this code available on research gate. I think it could help you.

Answer 2

[lu zhao]

I tweaked your code with RK and you see if that's correct.

[lu zhao]

% parameters of sinusoidal displacement and velocity

Ampl = 0.013;

freqHz = 0.6;

t=0:1e-3:10;

h=1e-3

% vectors describing virtual displacement and virtual velocity

X_virt = Ampl*sin(2*pi*freqHz*t);

V_virt = Ampl*2*pi*freqHz*cos(2*pi*freqHz*t);

% vector of constant parameters

vp = [650 300 8e4 0 1.2e-6 1e-6 15 2];

% solving differential equation

% f=@(t,z)vp(7).*V_virt-vp(5).*abs(V_virt).*z.*(abs(z)).^(vp(8)-1)-vp(6).*V_virt.*(abs(z)).^(vp(8))

z=0

n=length(t)

for i=1:n-1

f=@(t,z)vp(7).*V_virt(i)-vp(5).*abs(V_virt(i)).*z.*(abs(z)).^(vp(8)-1)-vp(6).*V_virt(i).*(abs(z)).^(vp(8))

k1=f(t(i),z(i))

k2=f(t(i)+0.5*h,z(i)+0.5*h*k1)

k3=f(t(i)+0.5*h,z(i)+0.5*h*k2)

k4=f(t(i)+h,z(i)+h*k3)

z(i+1)=z(i)+h/6.*(k1+2.*k2+2.*k3+k4)

end

Fmodel = vp(1)*V_virt + vp(2)*X_virt + vp(3)*z ;

% dx = V_virt;

% syms z(dx)

% z(dx) = dsolve (diff(z) == -vp(5)*abs(dx)*z*(abs(z))^(vp(8)-1)...

% -vp(6)*dx*(abs(z))^(vp(8)) + vp(7)*dx)

% % computing FModel

% Fmodel = vp(1)*V_virt + vp(2)*X_virt + vp(3)*z(dx) + vp(4);

% plotting

figure('Name', 'Test BW');

grid on;

plot (V_virt, Fmodel)