what is wrong in my code , there is a mistake in it because it is not converged
1 view (last 30 days)
Show older comments
clear;clc;warning off;
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta wx w
Kx=1500;Ky=1500;Kxt=150;Kyt=150;m=10;mt=1;
Lx=1.;
Ly=1.;
wx=(Kx/m)^0.5;
%mt is for TMD
Cx=2*0.02*(Kx/m)^0.5*m; Cy=2*0.02*(Ky/m)^0.5*m; Cxt=2*0.02*(Kxt/mt)^0.5*mt; Cyt=2*0.02*(Kyt/mt)^0.5*mt;
Beta=pi/6;
w=1;
%--------------------------------------------------------------------------
A=.03;
dt=.01;tf=30.;t=0:dt:tf;n=tf/dt;%tsp=time step; tf=final time;
load('CHICHI0968.mat');%CHICHI0968.mat-max of elcentro=0.3487
%ug=9.81*CHICHI0968(2,1:n+1);%%ELCENTRO_NS0348;A*(wx*w)^2*sin(wx*w*t);
ug=A*(wx*w)^2*sin(wx*w*t);
xg=ug;
%--------------------------------------------------------------------------
x0N=[0 0 0 0 0 0 0 0];x0L=x0N;xjN(1,:)=x0N;xjL(1,:)=x0L;
for j=1:n
%for j=1:3500
tint=dt*[j-1 j];%tint=time interval
[tN,xN] = ode45(@nonlinearmodel,tint,x0N);
xjN(j+1,:)=xN(length(tN),:);
x0N=xN(length(tN),:);
[tL,xL] = ode45(@linearmodel,tint,x0L);
xjL(j+1,:)=xL(length(tL),:);
x0L=xL(length(tL),:);
j;
end
%--------------------------------------------------------------------------
ux1=[xjN(:,1) xjL(:,1)];ux2=[xjN(:,2) xjL(:,2)];
uy1=[xjN(:,3) xjL(:,3)];uy2=[xjN(:,4) xjL(:,4)];
%ux1=[xjL(:,1)];ux2=[xjL(:,2)];
%uy1=[xjN(:,3) xjL(:,3)];uy2=[xjN(:,4) xjL(:,4)];
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta wx
% Define parameters
Kx=1500; Ky=1500; Kxt=150; Kyt=150; m=10; mt=1;
Lx=1.; Ly=1.;
wx=(Kx/m)^0.5;
Cx=2*0.02*(Kx/m)^0.5*m; Cy=2*0.02*(Ky/m)^0.5*m; Cxt=2*0.02*(Kxt/mt)^0.5*mt; Cyt=2*0.02*(Kyt/mt)^0.5*mt;
Beta=pi/6;
w=1;
% Define optimization problem
fun = @(A) cost_function(A);
x0 = 0.03;
lb = 0.01;
ub = 0.1;
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[A_opt, J_opt] = fmincon(fun, x0, [], [], [], [], lb, ub, [], options);
figure;
subplot(2,1,1);
plot(t, ux1(:,1), 'b', t, ux1(:,2), 'r');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Nonlinear Model', 'Linear Model');
title('Displacement of Main Mass (ux1)');
subplot(2,1,2);
plot(t, uy1(:,1), 'b', t, uy1(:,2), 'r');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Nonlinear Model', 'Linear Model');
title('Displacement of TMD (uy1)');
% Define function for cost function
function J = cost_function(A)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta wx w
dt=0.01;
tf=30.;
t=0:dt:tf;
n=tf/dt;
load('CHICHI0968.mat');
ug=A*(wx*w)^2*sin(wx*w*t);
xg=ug;
x0N=[0 0 0 0 0 0 0 0];
x0L=x0N;
xjN(1,:)=x0N;
xjL(1,:)=x0L;
for j=1:n
tint=dt*[j-1 j];
[tN,xN] = ode45(@nonlinearmodel,tint,x0N);
xjN(j+1,:)=xN(length(tN),:);
x0N=xN(length(tN),:);
[tL,xL] = ode45(@linearmodel,tint,x0L);
xjL(j+1,:)=xL(length(tL),:);
x0L=xL(length(tL),:);
end
ux1=[xjN(:,1) xjL(:,1)];
ux2=[xjN(:,2) xjL(:,2)];
uy1=[xjN(:,3) xjL(:,3)];
uy2=[xjN(:,4) xjL(:,4)];
% Define cost function as the maximum displacement of the main structure
J = max(abs(ux1(:,1)));
end
% Define differential equation for nonlinear model
function dx=nonlinearmodel(t,x)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta
dx = zeros(8,1);
dx(1)=x(5);
dx(2)=x(6);
dx(3)=x(7);
dx(4)=x(8);
dx(5)=1/m*(-Cx*x(5)-Kx*x(1)-Cxt*(x(5)-x(6))+Kxt*((x(2)-x(1))+(x(4)-x(3))^2/(2*Lx)-(x(2)-x(1))*(x(4)-x(3))^2/Lx^2))...
+1/m *Kyt*((x(2)-x(1))*(x(4)-x(3))/Ly -(x(2)-x(1))*(x(4)-x(3))^2/Ly^2+(x(2)-x(1))^3/(2*Ly^2))-xg(j)*cos(Beta);
dx(6)=1/mt*(Cxt*(x(5)-x(6))-Kxt*((x(2)-x(1))+(x(4)-x(3))^2/(2*Lx)-(x(2)-x(1))*(x(4)-x(3))^2/Lx^2))-xg(j)*cos(Beta);
dx(7)=1/m*(-Cy*x(7)-Ky*x(3)-Cyt*(x(7)-x(8))+Kyt*((x(4)-x(3))+(x(2)-x(1))^2/(2*Ly)-(x(4)-x(3))*(x(2)-x(1))^2/Ly^2))...
+1/m*Kxt*((x(4)-x(3))*(x(2)-x(1))/Lx -(x(4)-x(3))*(x(2)-x(1))^2/Lx^2+(x(4)-x(3))^3/(2*Lx^2))-xg(j)*sin(Beta);
dx(8)=1/mt*(Cxt*(x(7)-x(8))-Kyt*((x(4)-x(3))+(x(2)-x(1))^2/(2*Ly)-(x(4)-x(3))*(x(2)-x(1))^2/Ly^2))...
-1/mt*Kxt*((x(4)-x(3))*(x(2)-x(1))/Lx -(x(4)-x(3))*(x(2)-x(1))^2/Lx^2+(x(4)-x(3))^3/(2*Lx^2))-xg(j)*sin(Beta);
end
% Define differential equation for linear model
function dx=linearmodel(t,x)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta
dx = zeros(8,1);
dx(1)=x(5);
dx(2)=x(6);
dx(3)=x(7);
dx(4)=x(8);
dx(5)=1/m*(-Cx*x(5)-Kx*x(1)-Cxt*(x(5)-x(6))+Kxt*((x(2)-x(1))))-xg(j)*cos(Beta);
dx(6)=1/mt*(Cxt*(x(5)-x(6))-Kxt*((x(2)-x(1))))-xg(j)*cos(Beta);
dx(7)=1/m*(-Cy*x(7)-Ky*x(3)-Cyt*(x(7)-x(8))+Kyt*((x(4)-x(3))))-xg(j)*sin(Beta);
dx(8)=1/mt*(Cxt*(x(7)-x(8))-Kxt*(x(4)-x(3)))-xg(j)*sin(Beta);
end
2 Comments
Answers (0)
See Also
Categories
Find more on Nonlinear Control in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!