Clear Filters
Clear Filters

i want some help in my code to do the optimization to the main structure and how to derive the max. response under harmonic excitation with different freq. and optimize it

2 views (last 30 days)
Doha Ali
Doha Ali on 26 May 2023
Edited: Doha Ali on 26 May 2023
clear;clc; warning off;
global Kx Ky Kz m Cx Cy Cz Kxt Kyt Kzt mt Cxt Cyt Czt xg L1x L2x L1y L2y L1z L2z j Beta
wx=2*pi/1.11;
wy=wx;
wz=wx*100;
m=100;
mt=5/100;
Kx=m*wx^2;
Ky=m*wy^2;
Kz=m*wz^2;
Cx=2*0.05*m*wx;
Cy=2*0.05*m*wy;
Cz=2*0.05*m*wz;
sigmax=0.01;
sigmay=0.001;
wxt=4.1/2;
wxt=wx;
wyt=wxt;
wzt=2.9*wxt+sigmay;
%wzt=32^0.5*wxt;
%wxt=3;
%wyt=wxt;
%wyt=4.5;
%mt=0.6;
whar=wxt;
%whar=wyt;
Kxt=mt*wxt^2;
%Kxt=21.6;
%Kyt=21.6/4*2.25;
Kyt=mt*wyt^2;
Kzt=mt*wzt^2;
L1x=0.17;
L2x=0.375-L1x;
L1y=0.17;
L2y=0.375-L1x;
L1z=0.215;
L2z=mt*9.81/Kzt;
%L2z=0.001;
L2=0.01;
L1=1.0;
%wxt=6;
%Lx=1.;
%Lz=1. -mt*9.81/Kyt;
Beta=0.0;
%mt is for TMD
Kisay=0.05;
%Kisay=0.020;
%mt is for TMD
Cxt=2*0.023*(Kxt/mt)^0.5*mt; Cyt=2*0.023*(Kyt/mt)^0.5*mt;Czt=2*0.018*(Kzt/mt)^0.5*mt;
w=1.39;
%--------------------------------------------------------------------------
A=2.1*0.2;
dt=.02;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*sin(wxt*w*t);
xg=ug;
%--------------------------------------------------------------------------
x0N=[0.014 0 0 0 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_structure,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)];ux2=[xjN(:,2)];ux3=[xjN(:,3)];
save('Nolin.txt','xjN','-ASCII');
%ux1=[xjL(:,1)];ux2=[xjL(:,2)];
%uy1=[xjN(:,3) xjL(:,3)];uy2=[xjN(:,4) xjL(:,4)];
% clear;clc;
global Kx Ky Kz m Cx Cy Cz Kxt Kyt Kzt mt Cxt Cyt Czt xg L1x L2x L1y L2y L1z L2z j Beta
% Define optimization parameters
nvars = 6; % Number of variables to be optimized
lb = [0.01 0.01 0.01 0.01 , 0.01]; % Lower bounds for the optimization variables
ub = [10 10 10 10 10 10 ]; % Upper bounds for the optimization variables
% Define the GA optimization options
options = optimoptions('ga', 'Display', 'iter', 'PlotFcn' ...
, @gaplotbestf, 'MaxGenerations', 6, 'PopulationSize', 10);
% Define the objective function to be minimized
objFunc = @(x) simulate_structure(x);
% Run the GA optimization
[x_opt, fval_opt] = ga(objFunc, nvars, [], [], [], [], lb, ub, [], options);
% Display the optimized parameters and the corresponding objective function value
disp('Optimized parameters:');
disp(x_opt);
disp('Objective function value:');
disp(fval_opt);
filename = 'Nolin.txt';
formatSpec = '%16f%16f%16f%16f%16f%16f%16f%16f%16f%16f%16f%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'TextType', 'string', 'ReturnOnError', false);
fclose(fileID);
Nolin = table(dataArray{1:end-1}, 'VariableNames', {'e02','e00','e1','e2','e3','e4','e5','e6','e7','e8','e9','e10'});
clearvars filename formatSpec fileID dataArray ans;
figure
dt = 0.02;
tf = 30.;
t = 0:dt:tf;
subplot(2, 1, 1)
plot(t, Nolin.e2);
subplot(2, 1, 2)
plot(t, x_opt(1,4)*sin(4*w*t/0.6905)/100, t, Nolin.e2);
legend("optimized Model", "Non-linear Model")
function f = simulate_structure(x)
global Kx Ky Kz m Cx Cy Cz Kxt Kyt Kzt mt Cxt Cyt Czt xg L1x L2x L1y L2y L1z L2z j Beta
% Define the system parameters
wx = 2*pi/1.11;
wy = wx;
wz = wx*100;
m = 100;
mt = 5/100;
Kx = m*wx^2;
Ky = m*wy^2;
Kz = m*wz^2;
Cx = 2*0.05*m*wx;
Cy = 2*0.05*m*wy;
Cz = 2*0.05*m*wz;
sigmax = 0.01;
sigmay = 0.001;
wxt = 4.1/2;
wxt = wx;
wyt = wxt;
wzt = 2.9*wxt+sigmay;
Kyt = mt*wyt^2;
Kzt = mt*wzt^2;
L1x = 0.17;
L2x = 0.375-L1x;
L1y = 0.17;
L2y = 0.375-L1x;
L1z = 0.215;
L2z = mt*9.81/Kzt;
L2 = 0.01;
L1 = 1.0;
Beta = 0;
Kisay = 0.05;
Cxt = 2*0.023*(Kxt/mt)^0.5*mt;
Cyt = 2*0.023*(Kyt/mt)^0.5*mt;
Czt = 2*0.018*(Kzt/mt)^0.5*mt;
w = 1.39;
A = 2.1*0.2;
dt = 0.1;
tf = 30.;
t = 0:dt:tf;
n = tf/dt;
ug = A*sin(wxt*w*t);
xg = ug;
% Update the system parameters with the optimized values
Kxt = x(1);
Kyt = x(2);
Kzt = x(3);
L1x = x(4);
L1y = x(5);
% Initialize the simulation
x0N = [0.014 0 0 0 0 0 0 0 0 0 0 0];
x0N(1,:) = x0N;
for j = 1:n
tint = dt*[j-1 j];
[tN,xN] = ode45(@nonlinearmodel_structure,tint,x0N);
xjN(j+1,:) = xN(length(tN),:);
x0N = xN(length(tN),:);
end
% Calculate the objective function value as the maximum displacement of the structure
% from its initial position over the simulation time
f = max(xjN(:,1));
end
function dx=nonlinearmodel_structure(t,x)
global Kx Ky Kz m Cx Cy Cz Kxt Kyt Kzt mt Cxt Cyt Czt xg L1x L2x L1y L2y L1z L2z j Beta
dx = zeros(12,1);
dx(1)=x(7);
dx(2)=x(8);
dx(3)=x(9);
dx(4)=x(10);
dx(5)=x(11);
dx(6)=x(12);
%dx(1)=x(4);
%dx(2)=x(5);
%dx(3)=x(6);
%Vatar_x=((Lx+xx)^2+yy^2)^0.5;
% Vatar_y=((yy+mt*9.81/Kyt+Ly)^2+xx^2)^0.5;
dx(7)=1/m*(-Cx*x(7)-Kx*x(1)+...
-(-Cxt*(x(10)-x(7))...
-Kxt/2*( (x(4)-x(1))+L2x+(x(6)-x(3))^2*0.5*L1x/(L1x+L2x)^2+(x(5)-x(2))^2*0.5*L1x/(L1x+L2x)^2-(x(4)-x(1))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(4)-x(1))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3)...
+Kxt/2*(-(x(4)-x(1))+L2x+(x(6)-x(3))^2*0.5*L1x/(L1x+L2x)^2+(x(5)-x(2))^2*0.5*L1x/(L1x+L2x)^2+(x(4)-x(1))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3+(x(4)-x(1))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3)...
-Kyt/2*((x(4)-x(1))* (x(5)-x(2)) *L1y/(L1y+L2y)^2+(x(4)-x(1))*L2y/(L1y+L2y)+0.5*(x(4)-x(1))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3- (x(4)-x(1)) *(x(5)-x(2))^2*L1y/(L1y+L2y)^3+0.5*(x(4)-x(1))^3*L1y/(L1y+L2y)^3)...
-Kyt/2*((x(4)-x(1))*(-(x(5)-x(2)))*L1y/(L1y+L2y)^2+(x(4)-x(1))*L2y/(L1y+L2y)+0.5*(x(4)-x(1))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3-(x(4)-x(1))*(-(x(5)-x(2)))^2*L1y/(L1y+L2y)^3+0.5*(x(4)-x(1))^3*L1y/(L1y+L2y)^3)...
-Kzt*((x(4)-x(1))*(x(6)-x(3))*L1z/(L1z+L2z)^2+(x(4)-x(1))*L2z/(L1z+L2z)+0.5*(x(4)-x(1))*(x(5)-x(2))^2*L1z/(L1z+L2z)^3-(x(4)-x(1))*(x(6)-x(3))^2*L1z/(L1z+L2z)^3+0.5*(x(4)-x(1))^3*L1z/(L1z+L2z)^3)))...
-xg(j)*cos(Beta);
dx(8)=1/m*(-Cy*x(8)-Ky*x(2)+...
-(-Cyt*(x(11)-x(8))...
-Kyt/2*( (x(5)-x(2))+L2y+(x(6)-x(3))^2*0.5*L1y/(L1y+L2y)^2+(x(4)-x(1))^2*0.5*L1y/(L1y+L2y)^2-(x(5)-x(2))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3-(x(5)-x(2))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3)...
+Kyt/2*(-(x(5)-x(2))+L2y+(x(6)-x(3))^2*0.5*L1y/(L1y+L2y)^2+(x(4)-x(1))^2*0.5*L1y/(L1y+L2y)^2+(x(5)-x(2))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3+(x(5)-x(2))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3)...
-Kxt/2*((x(5)-x(2))* (x(4)-x(1))*L1x/(L1x+L2x)^2+(x(5)-x(2))*L2x/(L1x+L2x)+0.5*(x(5)-x(2))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(5)-x(2))* (x(4)-x(1))^2*L1x/(L1x+L2x)^3+0.5*(x(5)-x(2))^3*L1x/(L1x+L2x)^3)...
-Kxt/2*((x(5)-x(2))*(-(x(4)-x(1)))*L1x/(L1x+L2x)^2+(x(5)-x(2))*L2x/(L1x+L2x)+0.5*(x(5)-x(2))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(5)-x(2))*(-(x(4)-x(1)))^2*L1x/(L1x+L2x)^3+0.5*(x(5)-x(2))^3*L1x/(L1x+L2x)^3)...
-Kzt*((x(5)-x(2))*(x(6)-x(3))*L1z/(L1z+L2z)^2+(x(5)-x(2))*L2z/(L1z+L2z)+0.5*(x(5)-x(2))*(x(4)-x(1))^2*L1z/(L1z+L2z)^3-(x(5)-x(2))*(x(6)-x(3))^2*L1z/(L1z+L2z)^3+0.5*(x(5)-x(2))^3*L1z/(L1z+L2z)^3)))...
-xg(j)*sin(Beta);
dx(9)=1/m*(-Cz*x(9)-Kz*x(3)+...
-(-Czt*(x(12)-x(9))...
-Kzt*((x(6)-x(3))+L2z+(x(4)-x(1))^2*0.5*L1z/(L1z+L2z)^2+(x(5)-x(2))^2*0.5*L1z/(L1z+L2z)^2-(x(6)-x(3))*(x(4)-x(1))^2*L1z/(L1z+L2z)^3-(x(6)-x(3))*(x(5)-x(2))^2*L1z/(L1z+L2z)^3)...
-Kyt/2*((x(6)-x(3))* (x(5)-x(2))*L1y/(L1y+L2y)^2+(x(6)-x(3))*L2y/(L1y+L2y)+0.5*(x(6)-x(3))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3-(x(6)-x(3))* (x(5)-x(2))^2*L1y/(L1y+L2y)^3+0.5*(x(6)-x(3))^3*L1y/(L1y+L2y)^3)...
-Kyt/2*((x(6)-x(3))*(-(x(5)-x(2)))*L1y/(L1y+L2y)^2+(x(6)-x(3))*L2y/(L1y+L2y)+0.5*(x(6)-x(3))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3-(x(6)-x(3))*(-(x(5)-x(2)))^2*L1y/(L1y+L2y)^3+0.5*(x(6)-x(3))^3*L1y/(L1y+L2y)^3)...
-Kxt/2*((x(6)-x(3))* (x(4)-x(1))*L1x/(L1x+L2x)^2+(x(6)-x(3))*L2x/(L1x+L2x)+0.5*(x(6)-x(3))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3-(x(6)-x(3)) *(x(4)-x(1))^2*L1x/(L1x+L2x)^3+0.5*(x(6)-x(3))^3*L1x/(L1x+L2x)^3)...
-Kxt/2*((x(6)-x(3))*(-(x(4)-x(1)))*L1x/(L1x+L2x)^2+(x(6)-x(3))*L2x/(L1x+L2x)+0.5*(x(6)-x(3))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3-(x(6)-x(3))*(-(x(4)-x(1)))^2*L1x/(L1x+L2x)^3+0.5*(x(6)-x(3))^3*L1x/(L1x+L2x)^3)))...
-xg(j)*sin(Beta)*0.+mt/mt*9.81;
dx(10)=1/mt*(-Cxt*(x(10)-x(7))...
-Kxt/2*( (x(4)-x(1))+L2x+(x(6)-x(3))^2*0.5*L1x/(L1x+L2x)^2+(x(5)-x(2))^2*0.5*L1x/(L1x+L2x)^2-(x(4)-x(1))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(4)-x(1))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3)...
+Kxt/2*(-(x(4)-x(1))+L2x+(x(6)-x(3))^2*0.5*L1x/(L1x+L2x)^2+(x(5)-x(2))^2*0.5*L1x/(L1x+L2x)^2+(x(4)-x(1))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3+(x(4)-x(1))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3)...
-Kyt/2*((x(4)-x(1))* (x(5)-x(2)) *L1y/(L1y+L2y)^2+(x(4)-x(1))*L2y/(L1y+L2y)+0.5*(x(4)-x(1))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3- (x(4)-x(1)) *(x(5)-x(2))^2*L1y/(L1y+L2y)^3+0.5*(x(4)-x(1))^3*L1y/(L1y+L2y)^3)...
-Kyt/2*((x(4)-x(1))*(-(x(5)-x(2)))*L1y/(L1y+L2y)^2+(x(4)-x(1))*L2y/(L1y+L2y)+0.5*(x(4)-x(1))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3-(x(4)-x(1))*(-(x(5)-x(2)))^2*L1y/(L1y+L2y)^3+0.5*(x(4)-x(1))^3*L1y/(L1y+L2y)^3)...
-Kzt*((x(4)-x(1))*(x(6)-x(3))*L1z/(L1z+L2z)^2+(x(4)-x(1))*L2z/(L1z+L2z)+0.5*(x(4)-x(1))*(x(5)-x(2))^2*L1z/(L1z+L2z)^3-(x(4)-x(1))*(x(6)-x(3))^2*L1z/(L1z+L2z)^3+0.5*(x(4)-x(1))^3*L1z/(L1z+L2z)^3))...
-xg(j)*cos(Beta);
dx(11)=1/mt*(-Cyt*(x(11)-x(8))...
-Kyt/2*( (x(5)-x(2))+L2y+(x(6)-x(3))^2*0.5*L1y/(L1y+L2y)^2+(x(4)-x(1))^2*0.5*L1y/(L1y+L2y)^2-(x(5)-x(2))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3-(x(5)-x(2))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3)...
+Kyt/2*(-(x(5)-x(2))+L2y+(x(6)-x(3))^2*0.5*L1y/(L1y+L2y)^2+(x(4)-x(1))^2*0.5*L1y/(L1y+L2y)^2+(x(5)-x(2))*(x(6)-x(3))^2*L1y/(L1y+L2y)^3+(x(5)-x(2))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3)...
-Kxt/2*((x(5)-x(2))* (x(4)-x(1))*L1x/(L1x+L2x)^2+(x(5)-x(2))*L2x/(L1x+L2x)+0.5*(x(5)-x(2))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(5)-x(2))* (x(4)-x(1))^2*L1x/(L1x+L2x)^3+0.5*(x(5)-x(2))^3*L1x/(L1x+L2x)^3)...
-Kxt/2*((x(5)-x(2))*(-(x(4)-x(1)))*L1x/(L1x+L2x)^2+(x(5)-x(2))*L2x/(L1x+L2x)+0.5*(x(5)-x(2))*(x(6)-x(3))^2*L1x/(L1x+L2x)^3-(x(5)-x(2))*(-(x(4)-x(1)))^2*L1x/(L1x+L2x)^3+0.5*(x(5)-x(2))^3*L1x/(L1x+L2x)^3)...
-Kzt*((x(5)-x(2))*(x(6)-x(3))*L1z/(L1z+L2z)^2+(x(5)-x(2))*L2z/(L1z+L2z)+0.5*(x(5)-x(2))*(x(4)-x(1))^2*L1z/(L1z+L2z)^3-(x(5)-x(2))*(x(6)-x(3))^2*L1z/(L1z+L2z)^3+0.5*(x(5)-x(2))^3*L1z/(L1z+L2z)^3))...
-xg(j)*sin(Beta);
dx(12)=1/mt*(-Czt*(x(12)-x(9))...
-Kzt*((x(6)-x(3))+L2z+(x(4)-x(1))^2*0.5*L1z/(L1z+L2z)^2+(x(5)-x(2))^2*0.5*L1z/(L1z+L2z)^2-(x(6)-x(3))*(x(4)-x(1))^2*L1z/(L1z+L2z)^3-(x(6)-x(3))*(x(5)-x(2))^2*L1z/(L1z+L2z)^3)...
-Kyt/2*((x(6)-x(3))* (x(5)-x(2))*L1y/(L1y+L2y)^2+(x(6)-x(3))*L2y/(L1y+L2y)+0.5*(x(6)-x(3))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3-(x(6)-x(3))* (x(5)-x(2))^2*L1y/(L1y+L2y)^3+0.5*(x(6)-x(3))^3*L1y/(L1y+L2y)^3)...
-Kyt/2*((x(6)-x(3))*(-(x(5)-x(2)))*L1y/(L1y+L2y)^2+(x(6)-x(3))*L2y/(L1y+L2y)+0.5*(x(6)-x(3))*(x(4)-x(1))^2*L1y/(L1y+L2y)^3-(x(6)-x(3))*(-(x(5)-x(2)))^2*L1y/(L1y+L2y)^3+0.5*(x(6)-x(3))^3*L1y/(L1y+L2y)^3)...
-Kxt/2*((x(6)-x(3))* (x(4)-x(1))*L1x/(L1x+L2x)^2+(x(6)-x(3))*L2x/(L1x+L2x)+0.5*(x(6)-x(3))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3-(x(6)-x(3)) *(x(4)-x(1))^2*L1x/(L1x+L2x)^3+0.5*(x(6)-x(3))^3*L1x/(L1x+L2x)^3)...
-Kxt/2*((x(6)-x(3))*(-(x(4)-x(1)))*L1x/(L1x+L2x)^2+(x(6)-x(3))*L2x/(L1x+L2x)+0.5*(x(6)-x(3))*(x(5)-x(2))^2*L1x/(L1x+L2x)^3-(x(6)-x(3))*(-(x(4)-x(1)))^2*L1x/(L1x+L2x)^3+0.5*(x(6)-x(3))^3*L1x/(L1x+L2x)^3))...
-xg(j)*sin(Beta)*0.+mt/mt*9.81;
end
  2 Comments
Star Strider
Star Strider on 26 May 2023
Ran in:
Please go into some detail as to what you are doing.
If you are optimising something, what is it?
If you are fitting data, which one of your differential equations are you fitting it to?
LD = load('CHICHI0968.mat');
Data = LD.CHICHI0968;
t = Data(1,:);
s = Data(2,:);
figure
plot(t, s)
grid
xlabel('Time')
ylabel('Signal')
It is best not to use global variables. Instead use the approach described in Passing Extra Parameters.
.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Assembly in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!