Escape de oscilaciones Ensayo-Simulación

Ecuación de movimiento

Generales

g; (*Aceleración de la gravedad en m/(s^2)*)
n=5; (*Número de brazos considerados*)
mmag=5;(*Masa del imán de control*)
me(*=45.5+4*mmag*); (*Masa de la estructura y los imanes de control*)
m={m1,m2,m3,m4,me}; (*Lista de masas en orden antihorario. me corresponde a la estructura*)
rmas=8; (*Radio de las masas pendulares*)
rmag=15/2; (*Radio del imán*)
l={l1,l2,l3,l4,0};(*Lista de longitudes de brazos en orden antihorario. El brazo del cuerpo estructural se considera 0*)
lmag=102; (*Longitud del brazo hasta el centro del imán en mm*)
ϕ=0; (*Ángulo entre el brazo principal y el del centro de masa*)

Los momentos de inercia están dados por:

Imass=(∑_(i=1)^4)((1/2)m[[i]](rmas^2)+m[[i]]((l[[i]])^2));(*Inercia de las masas pendulares como cuerpos rígidos*)
Imass
Is=(∑_(i=1)^4)((1/4)mmag (rmag^2)+mmag (lmag^2))+347125;(*Inercia de los imanes y estructura plástica*)
Is//N

Imass+Is

Determinación experimental del momento de inercia

Mediante la determinación de la aceleración angular del péndulo desde su posición de equilibrio cuando se le agrega un torque dado por τ=(Subscript[P, I])*(Subscript[r, I]), es posible determinar la inercia de éste como la relación

τ=Iα->I=(((Subscript[P, I])(Subscript[r, I]))/α); [N m (s^2)]->[Kg (m^2)]

El modelo 3D de Solidworks entrega un momento de inercia por elementos finitos de 981884.83 (gmm^2). Utilizando las fórmulas obtenidas mediante la ecuación de Lagrange y para los valores de longitudes y masas medidos, el momento de inercia teórico es de 1027210 (gmm^2), es decir que existe una discrepancia entre estos del 4,6%. En cuanto al momento de inercia obtenido mediante el ensayo de aceleración angular, considerando una masa (Subscript[m, I])=37,7 g y cuyo centro de masa está a 85 mm del eje de rotación, el ensayo otorga una aceleración de

Carga de datos de Ensayos

(*Ensayo N°1*)
IEnsayo1=Import["D:\\Onedrive\\Beca Doctoral CONICET\\Papers y Congresos\\2023 - Control de N-pendulo por motor de pulsos\\Ensayos momento de inercia\\Ensayo 1.txt","Table"];
IE1vel={};
Do[AppendTo[IE1vel,20*(IEnsayo1[[i+1]]-IEnsayo1[[i]])(π/180)],{i,1,Length@IEnsayo1-1}]
IE1acel={};
Do[AppendTo[IE1acel,20*(IE1vel[[i+1]]-IE1vel[[i]])],{i,1,Length@IE1vel-1}]
IE1acel
(*Ensayo N°2*)
IEnsayo2=Import["D:\\Onedrive\\Beca Doctoral CONICET\\Papers y Congresos\\2023 - Control de N-pendulo por motor de pulsos\\Ensayos momento de inercia\\Ensayo 2.txt","Table"];
IE2vel={};
Do[AppendTo[IE2vel,20*(IEnsayo2[[i+1]]-IEnsayo2[[i]])(π/180)],{i,1,Length@IEnsayo2-1}]
IE2acel={};
Do[AppendTo[IE2acel,20*(IE2vel[[i+1]]-IE2vel[[i]])],{i,1,Length@IE2vel-1}]
IE2acel
(*Ensayo N°3*)
IEnsayo3=Import["D:\\Onedrive\\Beca Doctoral CONICET\\Papers y Congresos\\2023 - Control de N-pendulo por motor de pulsos\\Ensayos momento de inercia\\Ensayo 3.txt","Table"];
IE3vel={};
Do[AppendTo[IE3vel,20*(IEnsayo3[[i+1]]-IEnsayo3[[i]])(π/180)],{i,1,Length@IEnsayo3-1}]
IE3acel={};
Do[AppendTo[IE3acel,20*(IE3vel[[i+1]]-IE3vel[[i]])],{i,1,Length@IE3vel-1}]
IE3acel

Seleccionando los valores de aceleración cercanos a θ=0 y promediando se obtiene la aceleración (Overscript[α, -]) en 1/(s^2)

αprom=(1/6)(IE1acel[[2,1]]+IE1acel[[3,1]]
	+IE2acel[[2,1]]+IE2acel[[3,1]]
	+IE3acel[[2,1]]+IE3acel[[3,1]])

Con este valor, el momento de inercia resultante vale

Iexperimental=((37.7*9810*85)/(αprom))

El valor obtenido experimentalmente es de 1690700 (gmm^2). Se infiere que la diferencia obtenida surge al no considerar ni en el modelo matemático ni en el de elementos finitos los momentos de inercia debidos al encoder, eje y elementos móviles del rodamiento. Al comparar el valor obtenido con un valor obtenido mediante prueba y error para forzar al sistema a converger con el ensayo de decremento libre N°4, la diferencia entre ambos es del 12,2% ((Subscript[I, pruebaforzado])=1506250).

Energía Potencial

Se considera la masa de la estructura como una masa pendular más, pero centrada en el eje de revolución.

x=Table[l[[i]]*Sin[θ[t]-ϕ+(i-1)((2π)/(n-1))],{i,1,n}];
y=Table[l[[i]]*Cos[θ[t]-ϕ+(i-1)((2π)/(n-1))],{i,1,n}];
h=Table[lc-y[[i]],{i,1,n}];
V=Collect[g(∑_(i=1)^n)m[[i]]h[[i]],{g,lc,Cos[θ[t]],Sin[θ[t]]}]

Las coordenadas del centro de masa están dadas por

Mp=(∑_(i=1)^n)m[[i]];
xc=(1/Mp)(∑_(i=1)^n)m[[i]] x[[i]]//FullSimplify
yc=(1/Mp)(∑_(i=1)^n)m[[i]] y[[i]]//FullSimplify

Mientras que el brazo y altura del centro de masa valen

lc=FullSimplify[(((xc^2)+(yc^2))^(1/2))]
hc=lc-yc//Simplify

Con lo que la energía potencial puede ser reescrita considerando solamente el centro de masas como

V=Assuming[{m1+m2+m3+m4+me>0,l2 m2==l4 m4,l1 m1-l3 m3>0},g Mp hc//FullSimplify](*Verificada!!! con las consideraciones*)

Energía Cinética

Para la determinación de la energía cinética, se separa el sistema entre las masas y la estructura. La energía cinética de la estructura estará dada por

Te=(1/2)me((((X'[t]))^2)+(((Y'[t]))^2))+(1/2)Iz ((θ'[t])^2)

donde X’(t) e Y’(t) son las velocidades del forzamiento externo, e “Iz” es el momento de inercia polar de la estructura referido al eje de rotación.

Para la determinación de la expresión de la energía cinética de las masas pendulares, es necesario determinar la velocidad relativa de cada una de las masas con respecto a la del centro de masa. Si se define a esta última como (Subscript[v, θc]) y se descompone en sus componentes (Subscript[v, xc]) & (Subscript[v, yc]), la expresión general de la energía cinética estará dada por:

virel={v1rel,v2rel,v3rel,v4rel};
Mi=(∑_(i=1)^(n-1))m[[i]];(*Masa pendular total*)
vxc=Assuming[m1+m2+m3+m4+me>0,θ'[t]*lc*Cos[θ[t]]//FullSimplify]
vyc=Assuming[m1+m2+m3+m4+me>0,θ'[t]*lc*Sin[θ[t]]//FullSimplify]
Tm=(1/2)Mi((((Y'[t]+vyc))^2)+(((X'[t]+vxc))^2))+(1/2)(∑_(i=1)^(n-1))m[[i]]((virel[[i]])^2)//FullSimplify

donde la velocidad relativa es igual a la tasa de cambio de θ(t) multiplicada por dos vectores constantes. Como se debe cumplir que (Subscript[x, c])'=(Subscript[y, c])θ' y (Subscript[y, c])'=-(Subscript[x, c])θ', y del mismo modo  (Subscript[x, i])'=(Subscript[y, i])θ' y (Subscript[y, i])'=-(Subscript[x, i])θ', se pueden proponer las componentes de velocidad de cada partícula como

virel=θ'[t]*Table[{y[[i]]-yc,-x[[i]]+xc},{i,1,n-1}]; virel//MatrixForm
((virel[[1,1]])^2)+((virel[[1,2]])^2)/.{l1->1,l2->1,l3->1,l4->1,m1->2,m2->1,m3->1,m4->1,me->1}//ExpandAll//Simplify (*La velocidad relativa es (Subscript[K, i])θ'(t)*)

vrel={};
Do[AppendTo[vrel,((virel[[i,1]])^2)+((virel[[i,2]])^2)//ExpandAll//Simplify],{i,1,n-1}]
krel=(vrel/((θ'[t])^2))(*Expresiones específicas de los (Subscript[K, i])*)

La energía cinética de las masas vale entonces

Tm=(1/2)Mi((((Y'[t]+vyc))^2)+(((X'[t]+vxc))^2))+(1/2)((θ'[t])^2)(∑_(i=1)^(n-1))m[[i]]krel[[i]]/.{krel[[1]]->k1,krel[[2]]->k2,krel[[3]]->k3,krel[[4]]->k4,((((((l1 m1-l3 m3))^2)+(((l2 m2-l4 m4))^2))^(1/2))/(m1+m2+m3+m4+me))->Lc}//FullSimplify

y la energía cinética total vale

Ttot=Te+Tm

Acción de control

La diferencia con el estudio normal es que aquí se agrega una condición para que la acción de control no inicie hasta haber alcanzado un estacionario oscilatorio estable. La simulación es mucho más larga en este caso.

Modelo Sigmoide

Debido a los problemas de discontinuidades finitas en la integración numérica, se opta por implementar el mismo modelo de control mediante funciones sigmoides que cuentan con la ventaja de ser continuas y suaves. La implementación se realiza a través de una función dependiente de la posición angular θ(t), y otra a de la velocidad angular θ'(t). Partiendo de la serie de Fourier que representa la fuerza de control, se deben aplicar las reglas de accionamiento de la función en base a una ley de control expresada como una función ν(θ,θ'). Para el desarrollo de esta ley, se experimentó con el péndulo como un péndulo magnético, accionando el pulso electromagnético para rangos incrementales de ángulos de acción, hasta que la energía cinética aplicada es suficiente para alcanzar la rotación del péndulo. Bajo esta condición y para mejorar la eficiencia, una vez alcanzada la rotación, el sistema debe desactivar la acción de control.

FFtIθaprox7=0.9765709035985618` +0.6272412207035813` Cos[18 θ[t]]-0.4164160227041833` Cos[36 θ[t]]-0.009217893227308383` Cos[54 θ[t]]-0.05376352781178362` Cos[72 θ[t]]-0.05042019504839935` Cos[90 θ[t]]-0.01103290569019784` Cos[108 θ[t]]-0.03101998837072826` Cos[126 θ[t]];
tac=lmag*FFtIθaprox7 (*Torque de la acción de control en gmm*)

el control de accionamiento basado en la posición se define como la suma de cuatro sigmoides. La velocidad de activación se define a partir de una pendiente R. Las posiciones inicial y final de accionamiento se definen de manera simétrica a ambos lados de la posición de reposo del péndulo, y los valores de dichas posiciones se obtuvo de manera experimental para valores crecientes alcanzar la capacidad del péndulo magnético de obtener rotaciones

θini=1(π/180);
θfin=9(π/180);
R=200; (*ramp*)
θacSig=-(1/(1+(E^(-R(Mod[θ[t],2π,-π]+θini)))))+(1/(1+(E^(-R(Mod[θ[t],2π,-π]+θfin)))))+(1/(1+(E^(-R(Mod[θ[t],2π,-π]-θini)))))-(1/(1+(E^(-R(Mod[θ[t],2π,-π]-θfin)))))
ploθac=Plot[θacSig,{θ[t],-1.4θfin,1.4θfin},AspectRatio->1,ImageSize->UpTo[300],Axes->False,Frame->True,LabelStyle->Directive[FontFamily->"Cambria",18],FrameStyle->Directive[FontFamily->"Cambria",18],PlotStyle->{Dashed,RGBColor["#2196F3"],Thickness[0.01]},PlotLabel->"Position rule",FrameTicks->{{{0,1},None},{{θini,-θfin,θfin},None}},FrameLabel->{θ,ν[θ]},Exclusions->None,GridLinesStyle->Directive[Gray, Dashed],GridLines->{ {0},{0}}];

Para la velocidad, se realizaron simulaciones con rangos de aplicación mayores hasta alcanzar la rotación, momento en el cual se desactiva la acción de control.

ωini=0;
ωfin=2π*2;
R=200; (*ramp*)
ωacSig=(1/(1+(E^(-R(Θ[t]+ωini)))))-(1/(1+(E^(-R(Θ[t]+ωfin)))))+(1/(1+(E^(-R(Θ[t]-ωini)))))-(1/(1+(E^(-R(Θ[t]-ωfin)))));
ploωac=Plot[ωacSig,{Θ[t],-1.4ωfin,1.4ωfin},AspectRatio->1,ImageSize->UpTo[300],Axes->False,Frame->True,LabelStyle->Directive[FontFamily->"Cambria",18],FrameStyle->Directive[FontFamily->"Cambria",18],PlotStyle->{Dashed,RGBColor["#2196F3"],Thickness[0.01]},PlotLabel->"Position rule",FrameTicks->{{{0,1},None},{{ωini,-ωfin,ωfin},None}},FrameLabel->{Θ,ν[Θ]},Exclusions->None,GridLinesStyle->Directive[Gray, Dashed],GridLines->{ {0},{0}}];

Para evitar las situaciones donde el accionamiento del dispositivo disminuya la energía total, se desarrolla una regla en la que el signo de la velocidad angular debe ser opuesto al de la posición angular para obtener un accionamiento. Esto es debido a que de no cumplirse, el mecanismo que desarrolla una fuerza de repulsión se acciona cuando la masa principal acaba de pasar por la posición de reposo, y el torque aplicado tiende a repeler hacia la posición de reposo.

condb=If[Sign[Mod[θ[t],2π,-π]]==Sign[Θ[t]],0,1];

La función del torque aplicado finalmente toma la forma

νac=FFtIθaprox7*θacSig*ωacSig*condb;
νacR=Piecewise[{{0,t<=700},{νac,t>700}}];

Ensayos

Ensayo1=Import["C:\\Users\\User\\H=45 w=1.25\\Osci2.txt","Table"];(*Osci-Osci*)
Ensayo2=Import["C:\\Users\\User\\H=45 w=1.5\\Osci1.txt","Table"];(*Osci-Rot*)

List1=Table[{Ensayo1[[i,1]]-Ensayo1[[1,1]],Ensayo1[[i,2]]*(π/180),((Ensayo1[[i+1,2]]-Ensayo1[[i,2]])/(Ensayo1[[i+1,1]]-Ensayo1[[i,1]]))*(π/180)},{i,1,Length@Ensayo1-1}];
List2=Table[{Ensayo2[[i,1]]-Ensayo2[[1,1]],Ensayo2[[i,2]]*(π/180),((Ensayo1[[i+1,2]]-Ensayo1[[i,2]])/(Ensayo1[[i+1,1]]-Ensayo1[[i,1]]))*(π/180)},{i,1,Length@Ensayo1-1}];

Osci-Osci

Experimental

List1Phase=Select[List1[[All,{2,3}]],#[[1]] <=List1[[Length@Ensayo1-1,1]]&];
List1PhaseINI=Table[List1Phase[[i]],{i,1,300}];
List1PhaseFIN=Table[List1Phase[[i]],{i,3100,3400}];
List1ΘDesp=Table[{Ensayo1[[i,1]]-Ensayo1[[3100,1]],Ensayo1[[i,2]]*(π/180),((Ensayo1[[i+1,2]]-Ensayo1[[i,2]])/(Ensayo1[[i+1,1]]-Ensayo1[[i,1]]))*(π/180)},{i,3100,3400}];

Osci1INIExp=ListPlot[List1[[All,{1,3}]],Joined->True,PlotRange->{{0,10},Automatic},PlotRange->{Automatic,{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True,FrameLabel->{t,Θ},FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,200,2}],None}},PlotStyle->{Thickness[0.0045],RGBColor["#D1C4E9"]},AspectRatio->5/7,ImageSize->Round[72*(8/(2 2.54))], Epilog->{Dashed,Line[{{0,0},{20,0}}]}];

Osci1FINExp=ListPlot[List1ΘDesp[[All,{1,3}]],Joined->True,PlotRange->{{0,10},Automatic},PlotRange->{Automatic,{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True,FrameLabel->{t,None},FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,200,2}],None}},PlotStyle->{Thickness[0.0045],RGBColor["#5E35B1"]},AspectRatio->5/7,ImageSize->Round[72*(8/(2 2.54))], Epilog->{Dashed,Line[{{0,0},{180,0}}]}];

Osci1Comb=Show[Osci1INIExp,Osci1FINExp];

OsciOsciExp=ListPlot[{List1PhaseINI,List1PhaseFIN},PlotRange->{{-π,π},{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True,FrameLabel->{θ,None},FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,-π,π,π/2}],None}},PlotStyle->{{PointSize[0.0075],RGBColor["#D1C4E9"]},{PointSize[0.005],RGBColor["#5E35B1"]}},AspectRatio->5/7,ImageSize->Round[72*(8/(2 2.54))],Epilog->{Dashed,RGBColor["#E0E0E0"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black,Text[Style["a",FontFamily->"Palatino Linotype",15],{((7π)/8),-200}]}];

GraphicsRow[{Osci1Comb,OsciOsciExp},ImageSize->Round[72*(16/2.54)],Spacings->0]

Osci-Rot

Experimental

List2Phase=Table[{Ensayo2[[i,2]]*(π/180),((Ensayo2[[i+1,2]]-Ensayo2[[i,2]])/(Ensayo2[[i+1,1]]-Ensayo2[[i,1]]))*(π/180)},{i,1,Length@Ensayo2-1}];
List2PhaseINI=Table[List2Phase[[i]],{i,1,220}];
List2PhaseTRANS=Table[List2Phase[[i]],{i,220,3045}];
List2PhaseFIN=Table[List2Phase[[i]],{i,3045,3185}];

OsciRotExp=ListPlot[{List2PhaseTRANS,List2PhaseINI,List2PhaseFIN},PlotRange->{{-π,π},{-5,15}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True,FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,-π,π,π/2}],None}},PlotStyle->{{PointSize[0.0045],RGBColor["#FFCC80"]},{PointSize[0.0075],RGBColor["#D1C4E9"]},{PointSize[0.0075],RGBColor["#D81B60"]}},AspectRatio->2/3,ImageSize->Round[72*(10/2.54)],Epilog->{{Dashing[{0.025,0.02}],Thickness[0.0025],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}]},Text[Style["a",FontFamily->"Palatino Linotype",15],{((7π)/8),-200}]}]

OsciEstE2=ListPlot[Select[List2[[All,{1,2}]],#[[1]]<=32.4 &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#D1C4E9"]}];
TransE2=ListPlot[Select[List2[[All,{1,2}]],32.35<=#[[1]]<143.9 &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#FFCC80"]}];
RotEstE2=ListPlot[Select[List2[[All,{1,2}]],143.89<=#[[1]] &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#D81B60"]}];
Show[{OsciEstE2,TransE2,RotEstE2},Ticks->{Table[i,{i,0,180,15}],Table[i,{i,-π,π,π/2}]},ImageSize->Round[72*(40/2.54)],AspectRatio->1/8]

Simulación

Osci-Osci

tamaño=6.5;

Experimental

List1Phase=Select[List1[[All,{2,3}]],#[[1]] <=List1[[Length@Ensayo1-1,1]]&];
List1PhaseINI=Table[List1Phase[[i]],{i,1,300}];
List1PhaseFIN=Table[List1Phase[[i]],{i,3100,3400}];
List1ΘDesp=Table[{Ensayo1[[i,1]]-Ensayo1[[3100,1]],Ensayo1[[i,2]]*(π/180),((Ensayo1[[i+1,2]]-Ensayo1[[i,2]])/(Ensayo1[[i+1,1]]-Ensayo1[[i,1]]))*(π/180)},{i,3100,3400}];

Osci1INIExp=ListPlot[List1[[All,{1,3}]],Joined->True,PlotRange->{{0,10},{-10,10}},PlotRange->{Automatic,{-10,10}},Axes->False,Frame->True(*,FrameLabel->{None,(OverDot[θ])}*),FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,200,2}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{Thickness[0.0045],RGBColor["#D1C4E9"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))], Epilog->{Dashed,RGBColor["#757575"],Line[{{0,0},{20,0}}],Black,Text[Style["a",FontFamily->"Palatino Linotype",13],{9.5,-8}]}];

Osci1FINExp=ListPlot[List1ΘDesp[[All,{1,3}]],Joined->True,PlotRange->{{0,10},Automatic},PlotRange->{Automatic,{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,11],Axes->False,Frame->True,PlotStyle->{Thickness[0.0045],RGBColor["#5E35B1"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))], Epilog->{Dashed,Line[{{0,0},{180,0}}]}];

Osci1EXPComb=Show[{Osci1INIExp,Osci1FINExp},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

OsciOsciExp=ListPlot[{List1PhaseINI,List1PhaseFIN},PlotRange->{{-π,π},{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True(*,FrameLabel->{None,None}*),FrameTicks->{{None,None},{Table[i,{i,-π,π,π/2}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{{PointSize[0.0075],RGBColor["#D1C4E9"]},{PointSize[0.005],RGBColor["#5E35B1"]}},Epilog->{Dashed,RGBColor["#757575"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black,Text[Style["b",FontFamily->"Palatino Linotype",13],{((7π)/8),-8}]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

GraphicsRow[{Osci1EXPComb,OsciOsciExp},ImageSize->Round[72*((2*tamaño)/2.54)],Spacings->0]

OsciOsciExpSalto=ListPlot[{List1PhaseINI,List1PhaseFIN},Joined->{False,True},PlotRange->{{-π/8,π/8},{6.25,7.75}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,15],Axes->False,Frame->True(*,FrameLabel->{None,None}*),FrameTicks->{{Table[i,{i,0,10,0.5}],None},{Table[i,{i,-π,π,π/16}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{{PointSize[0.025],RGBColor["#D1C4E9"]},{PointSize[0.025],Thickness[0.01],RGBColor["#5E35B1"]}},Epilog->{Dashed,RGBColor["#757575"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]]

Numérico

g=9810; (*mm/(s^2)*)
l1=85;l2=85;l3=65;l4=85; (*Longitudes en mm*)
m1=61;m2=43.9;m3=23.7;m4=43.9;me=66.7; (*Pesos en g*)
m={m1,m2,m3,m4,me};
mP=(∑_(i=1)^(n-1))m[[i]]; (*Masa de los péndulos*)
Iz=347125.27; (*Izz de Solidworks, momento de inercia polar de la estructura*)
Imasas=krel[[1]] m1+krel[[2]] m2+krel[[3]] m3+krel[[4]] m4; (*Momento de inercia polar debido a las masas pendulares. Junto con Icm resultan en el momento de inercia total de las masas*)
TrOpt=294265; (*Optimizado con puntaje 984*)
COpt=71214; (*Optimizado con puntaje 984*)
Icorr=125000; (*Optimizado con puntaje 984*)
Icm=(lc^2) (m1+m2+m3+m4); (*Momento de inercia del centro de masa*)
KFac= 3.05 (10^7);
tact=300;
Fact=Piecewise[{{0,t<tact},{1,t>=tact}}];;
Hey=45(*Semiamplitud de forzamiento vertical*);
Ωey=2π*1.25(*Frecuencia del forzamiento vertical*); 
fY[t]=Hey Cos[Ωey t];
tmax=800;

θ0=-0.646189406896826368`4.(*Posición inicial*)
Θ0=-5.738525390625`4.(*Velocidad inicial*)
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
SolMovθ=NDSolve[{
θ'[t]==Θ[t],
(Imasas+Iz+Icm+Icorr) Θ'[t]+COpt Θ[t]+TrOpt Sign[Θ[t]]+(g(l1 m1-l3 m3)+lc mP D[fY[t],{t,2}])Sin[θ[t]]-Fact KFac νac==0,
θ[0]==θ0,Θ[0]==Θ0},
{θ,Θ},{t,0,tmax},MaxSteps->Infinity,Method->"StiffnessSwitching"];
θsolu=θ[t]/.SolMovθ;
Θsolu=Θ[t]/.SolMovθ;

ΘsoluINI=Table[{i-290,(Θsolu/.t->i)[[1]]},{i,290,300,0.001}];
ΘsoluFIN=Table[{i-tmax+11,(Θsolu/.t->i)[[1]]},{i,tmax-11,tmax,0.001}];

GraphIterativoΘINI=ListLinePlot[ΘsoluINI,PlotRange->{{0,10},{-10,10}},Axes->False,Frame->True(*,FrameLabel->{t,(OverDot[θ])}*),RotateLabel->False,FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,800,2}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{Thickness[0.0045],RGBColor["#D1C4E9"]}, Epilog->{Dashed,RGBColor["#757575"],Line[{{0,0},{180,0}}],Black,Text[Style["c",FontFamily->"Palatino Linotype",13],{9.5,-8}]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

GraphIterativoΘFIN=ListLinePlot[ΘsoluFIN,PlotRange->{Automatic,{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,11],Axes->False,Frame->True,PlotStyle->{Thickness[0.0045],RGBColor["#5E35B1"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

Osci1NUMComb=Show[{GraphIterativoΘINI,GraphIterativoΘFIN},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];
(*Retatro de fase*)

listPhasePortraitINI={};
Do[AppendTo[listPhasePortraitINI,{Mod[θsolu,2π,-π][[1]]/.{t->i,θ[t]->Mod[θ[t],2π,-π]},Θsolu[[1]]/.t->i}],{i,290,300,0.005}]
listPhasePortraitFIN={};
Do[AppendTo[listPhasePortraitFIN,{Mod[θsolu,2π,-π][[1]]/.{t->i,θ[t]->Mod[θ[t],2π,-π]},Θsolu[[1]]/.t->i}],{i,tmax-30,tmax,0.005}]

GraphEF=ListPlot[{listPhasePortraitINI,listPhasePortraitFIN},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/2.54)],PlotRange->{{-π,π},{-10,10}},FrameTicks->{{None,None},{Table[i,{i,-π,π,π/2}],None}},Axes->False,Frame->True(*,FrameLabel->{None,None}*),LabelStyle->{FontFamily->"Palatino Linotype",11},RotateLabel->False,Joined->False,PlotStyle->{{PointSize[0.00025],Opacity[1],RGBColor["#D1C4E9"]},{PointSize[0.00025],Opacity[1],RGBColor["#5E35B1"]}},Epilog->{Dashed,RGBColor["#757575"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black,Text[Style["d",FontFamily->"Palatino Linotype",13],{((7π)/8),-8}]}];(*Hacer como lista y no unir los puntos*)

OsciOsciNum=GraphicsRow[{Osci1NUMComb,GraphEF},ImageSize->Round[72*((2*tamaño)/2.54)],Spacings->0]

OsciOsciNumSalto=ListPlot[{listPhasePortraitINI,listPhasePortraitFIN},Joined->{False,True},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/2.54)],PlotRange->{{-π/8,π/8},{6.25,7.75}},FrameTicks->{{Table[i,{i,0,10,0.5}],None},{Table[i,{i,-π,π,π/16}],None}},Axes->False,Frame->True(*,FrameLabel->{None,None}*),LabelStyle->{FontFamily->"Palatino Linotype",11},RotateLabel->False,Joined->False,PlotStyle->{{PointSize[0.025],Opacity[1],RGBColor["#D1C4E9"]},{PointSize[0.025],Thickness[0.01],Opacity[1],RGBColor["#5E35B1"]}},Epilog->{Dashed,RGBColor["#757575"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black}]

Combinación

OsciOsciComb=GraphicsGrid[{{Osci1EXPComb,OsciOsciExp},{Osci1NUMComb,GraphEF}},ImageMargins->{{5,0},{15,0}},Spacings->{-15,-5},
Epilog->{Text[Style["(OverDot[θ])",FontFamily->"Palatino Linotype",11],{0,-62.5}],Text[Style["(OverDot[θ])",FontFamily->"Palatino Linotype",11],{0,-205}],Text[Style["t",FontFamily->"Palatino Linotype",11],{90,-280}],Text[Style["θ",FontFamily->"Palatino Linotype",11],{242,-280}]}]
Export["OsciOsci.pdf",OsciOsciComb]

OsciOsciJump=GraphicsGrid[{{OsciOsciExpSalto,OsciOsciNumSalto}},
Epilog->{Text[Style["(OverDot[θ])",FontFamily->"Palatino Linotype",11],{9,-80}],Text[Style["θ",FontFamily->"Palatino Linotype",11],{300,-145}]}]
Export["OsciOsciJump.pdf",OsciOsciJump]

Osci-Rot

Experimental

List2Phase=Table[{Ensayo2[[i,2]]*(π/180),((Ensayo2[[i+1,2]]-Ensayo2[[i,2]])/(Ensayo2[[i+1,1]]-Ensayo2[[i,1]]))*(π/180)},{i,1,Length@Ensayo2-1}];
List2PhaseINI=Table[List2Phase[[i]],{i,1,220}];
List2PhaseTRANS=Table[List2Phase[[i]],{i,220,3045}];
List2PhaseFIN=Table[List2Phase[[i]],{i,3045,3185}];

List2ΘDesp=Table[{Ensayo2[[i,1]]-Ensayo2[[3100,1]],Ensayo2[[i,2]]*(π/180),If[Abs[((Ensayo2[[i+1,2]]-Ensayo2[[i,2]])/(Ensayo2[[i+1,1]]-Ensayo2[[i,1]]))*(π/180)]<20,((Ensayo2[[i+1,2]]-Ensayo2[[i,2]])/(Ensayo2[[i+1,1]]-Ensayo2[[i,1]]))*(π/180),((Ensayo2[[i,2]]-Ensayo2[[i-1,2]])/(Ensayo2[[i,1]]-Ensayo2[[i-1,1]]))*(π/180)]},{i,3100,3400}];

Osci2INIExp=ListPlot[List2[[All,{1,3}]],Joined->True,PlotRange->{{0,10},{-10,15}},Axes->False,Frame->True(*,FrameLabel->{None,(OverDot[θ])}*),FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,200,2}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{Thickness[0.0045],RGBColor["#D1C4E9"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))], Epilog->{Dashed,RGBColor["#757575"],Line[{{0,0},{20,0}}],Black,Text[Style["a",FontFamily->"Palatino Linotype",13],{9.5,-8}]}];

Osci2FINExp=ListPlot[List2ΘDesp[[All,{1,3}]],Joined->True,AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,11],Axes->False,Frame->True,PlotStyle->{Thickness[0.0045],RGBColor["#D81B60"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))], Epilog->{Dashed,Line[{{0,0},{180,0}}]}];
Osci2EXPComb=Show[{Osci2INIExp,Osci2FINExp},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

OsciRotExp=ListPlot[{List2PhaseTRANS,List2PhaseINI,List2PhaseFIN},PlotRange->{{-π,π},{-10,15}},FrameStyle->Directive[FontFamily->"Palatino Linotype",Black,11],Axes->False,Frame->True,FrameTicks->{{None,None},{Table[i,{i,-π,π,π/2}],None}},PlotStyle->{{PointSize[0.0045],RGBColor["#FFCC80"]},{PointSize[0.0075],RGBColor["#D1C4E9"]},{PointSize[0.0075],RGBColor["#D81B60"]}},AspectRatio->2/3,ImageSize->Round[72*((tamaño)/2.54)],Epilog->{{Dashing[{0.025,0.02}],Thickness[0.0025],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}]},Text[Style["b",FontFamily->"Palatino Linotype",13],{((7π)/8),-8}]}];

OsciRotEXP=GraphicsRow[{Osci2EXPComb,OsciRotExp},ImageSize->Round[72*((2*tamaño)/2.54)],Spacings->0]

OsciEstE2=ListPlot[Select[List2[[All,{1,2}]],#[[1]]<=32.4 &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#D1C4E9"]}];
TransE2=ListPlot[Select[List2[[All,{1,2}]],32.35<=#[[1]]<143.9 &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#FFCC80"]}];
RotEstE2=ListPlot[Select[List2[[All,{1,2}]],143.89<=#[[1]] &],Joined->True,PlotRange->{{0,165},{-π,π}},PlotStyle->{PointSize[0.0015],RGBColor["#D81B60"]}];
Show[{OsciEstE2,TransE2,RotEstE2},Ticks->{Table[i,{i,0,180,15}],Table[i,{i,-π,π,π/2}]},ImageSize->Round[72*(40/2.54)],AspectRatio->1/8]

Numérico

g=9810; (*mm/(s^2)*)
l1=85;l2=85;l3=65;l4=85; (*Longitudes en mm*)
m1=61;m2=43.9;m3=23.7;m4=43.9;me=66.7; (*Pesos en g*)
m={m1,m2,m3,m4,me};
mP=(∑_(i=1)^(n-1))m[[i]]; (*Masa de los péndulos*)
Iz=347125.27; (*Izz de Solidworks, momento de inercia polar de la estructura*)
Imasas=krel[[1]] m1+krel[[2]] m2+krel[[3]] m3+krel[[4]] m4; (*Momento de inercia polar debido a las masas pendulares. Junto con Icm resultan en el momento de inercia total de las masas*)
TrOpt=294265; (*Optimizado con puntaje 984*)
COpt=71214; (*Optimizado con puntaje 984*)
Icorr=800000; (*Optimizado con puntaje 984*)
Icm=(lc^2) (m1+m2+m3+m4); (*Momento de inercia del centro de masa*)
KFac= 3.05 (10^7);
tact=300;
Fact=Piecewise[{{0,t<tact},{1,t>=tact}}];;
Hey=45(*Semiamplitud de forzamiento vertical*);
Ωey=2π*1.25(*Frecuencia del forzamiento vertical*); 
fY[t]=Hey Cos[Ωey t];
tmax=2000;

θ0=-0.646189406896826368`4.(*Posición inicial*)
Θ0=-5.738525390625`4.(*Velocidad inicial*)
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
SolMovθ=NDSolve[{
θ'[t]==Θ[t],
(Imasas+Iz+Icm+Icorr) Θ'[t]+COpt Θ[t]+TrOpt Sign[Θ[t]]+(g(l1 m1-l3 m3)+lc mP D[fY[t],{t,2}])Sin[θ[t]]-Fact KFac νac==0,
θ[0]==θ0,Θ[0]==Θ0},
{θ,Θ},{t,0,tmax},MaxSteps->Infinity,Method->"StiffnessSwitching"];
θsolu=θ[t]/.SolMovθ;
Θsolu=Θ[t]/.SolMovθ;

ΘsoluINI=Table[{i-290,(Θsolu/.t->i)[[1]]},{i,290,300,0.001}];
ΘsoluFIN=Table[{i-tmax+11,(Θsolu/.t->i)[[1]]},{i,tmax-11,tmax,0.001}];

GraphIterativoΘINI=ListLinePlot[ΘsoluINI,PlotRange->{{0,10},{-10,15}},Axes->False,Frame->True(*,FrameLabel->{t,(OverDot[θ])}*),RotateLabel->False,FrameTicks->{{Table[i,{i,-20,20,5}],None},{Table[i,{i,0,800,2}],None}},LabelStyle->{FontFamily->"Palatino Linotype",11},PlotStyle->{Thickness[0.0045],RGBColor["#D1C4E9"]}, Epilog->{Dashed,RGBColor["#757575"],Line[{{0,0},{180,0}}],Black,Text[Style["c",FontFamily->"Palatino Linotype",13],{9.5,-8}]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

GraphIterativoΘFIN=ListLinePlot[ΘsoluFIN,PlotRange->{Automatic,{-10,10}},AxesStyle->Directive[FontFamily->"Palatino Linotype",Black,11],Axes->False,Frame->True,PlotStyle->{Thickness[0.0045],RGBColor["#D81B60"]},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];

Osci1NUMComb=Show[{GraphIterativoΘINI,GraphIterativoΘFIN},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/( 2.54))]];
(*Retatro de fase*)

listPhasePortraitINI={};
Do[AppendTo[listPhasePortraitINI,{Mod[θsolu,2π,-π][[1]]/.{t->i,θ[t]->Mod[θ[t],2π,-π]},Θsolu[[1]]/.t->i}],{i,290,300,0.005}]
listPhasePortraitTRANS={};
Do[AppendTo[listPhasePortraitTRANS,{Mod[θsolu,2π,-π][[1]]/.{t->i,θ[t]->Mod[θ[t],2π,-π]},Θsolu[[1]]/.t->i}],{i,301,tmax-30,0.1}]
listPhasePortraitFIN={};
Do[AppendTo[listPhasePortraitFIN,{Mod[θsolu,2π,-π][[1]]/.{t->i,θ[t]->Mod[θ[t],2π,-π]},Θsolu[[1]]/.t->i}],{i,tmax-30,tmax,0.005}]

GraphEF=ListPlot[{listPhasePortraitTRANS,listPhasePortraitINI,listPhasePortraitFIN},AspectRatio->5/7,ImageSize->Round[72*((tamaño)/2.54)],PlotRange->{{-π,π},{-10,15}},FrameTicks->{{None,None},{Table[i,{i,-π,π,π/2}],None}},Axes->False,Frame->True(*,FrameLabel->{None,None}*),LabelStyle->{FontFamily->"Palatino Linotype",11},RotateLabel->False,Joined->False,PlotStyle->{{PointSize[0.00025],Opacity[1],RGBColor["#FFCC80"]},{PointSize[0.00025],Opacity[1],RGBColor["#D1C4E9"]},{PointSize[0.00025],Opacity[1],RGBColor["#D81B60"]}},Epilog->{Dashed,RGBColor["#757575"],Line[{{0,-400},{0,800}}],Line[{{-2π,0},{2π,0}}],Black,Text[Style["d",FontFamily->"Palatino Linotype",13],{((7π)/8),-8}]}];(*Hacer como lista y no unir los puntos*)

OsciOsciNum=GraphicsRow[{Osci1NUMComb,GraphEF},ImageSize->Round[72*((2*tamaño)/2.54)],Spacings->0]

Combinación

OsciRotComb=GraphicsGrid[{{Osci2EXPComb,OsciRotExp},{Osci1NUMComb,GraphEF}},ImageMargins->{{5,0},{15,0}},Spacings->{-15,-5},
Epilog->{Text[Style["(OverDot[θ])",FontFamily->"Palatino Linotype",11],{0,-64}],Text[Style["(OverDot[θ])",FontFamily->"Palatino Linotype",11],{0,-200}],Text[Style["t",FontFamily->"Palatino Linotype",11],{92.5,-270}],Text[Style["θ",FontFamily->"Palatino Linotype",11],{242.5,-270}]}]
Export["OsciRot.pdf",OsciRotComb]