Determinación de Parámetros

Necesarios

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\\2025 - 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\\2025 - 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\\2025 - 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

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=12;
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

Plot[FFtIθaprox7,{θ[t],-12 (π/180),12(π/180)}];
νac=FFtIθaprox7*θacSig*ωacSig*condb;
τac=νac*(lmag/1000);

ContoursColors={RGBColor["#d81b60"],RGBColor["#e05993"],RGBColor["#e285bd"],RGBColor["#e2addd"],RGBColor["#e7d2f2"],RGBColor["#ffffff"],RGBColor["#ffffff"],RGBColor["#d2ceea"],RGBColor["#b0aada"],RGBColor["#8c88cb"],RGBColor["#6768bb"],RGBColor["#3949ab"]};
Borde=RGBColor["#d4af37"];


FiguraC=ContourPlot[τac,{θ[t],-0.3 ,0.3},{Θ[t],-1.25 ωfin,1.25 ωfin},PlotRange->{Full,Full,{-1.5,1.5}},Contours->11,ColorFunction->ContoursColors,ContourStyle->None,FrameLabel->{"θ","θ'"},FrameStyle->Directive[FontFamily->"Palatino Linotype",Black,11],
FrameTicks->{{Table[i,{i,-15,15,5}],None},{Table[i,{i,-20 (π/180),20 (π/180),5 (π/180)}],None}},
PlotLegends->BarLegend[{Automatic,{-1.5,1.5}},LegendLabel->Text[Style["(Subscript[F, t])",FontFamily->"Palatino Linotype",11]]],LabelStyle->Directive[FontFamily->"Palatino Linotype",Black,11],AspectRatio->1,ImageSize->Round[((72*8)/2.54)],MaxRecursion->6,Exclusions->None];
ColorBlend[z_]:=Blend[{{-(1.5/100),RGBColor["#3949AB"]},{0,White},{(1.5/100),RGBColor["#D81B60"]}},z];
FiguraD=DensityPlot[τac,{θ[t],-0.2 ,0.2},{Θ[t],-1.1 ωfin,1.1 ωfin},PlotRange->{Full,Full,{-1.5,1.5}},ColorFunction->(ColorBlend[#]&),ColorFunctionScaling->False,MaxRecursion->5,Exclusions->None,FrameLabel->{"θ","(OverDot[θ])"},FrameStyle->Directive[FontFamily->"Palatino Linotype",Black,11],
FrameTicks->{{Table[i,{i,-15,15,5}],None},{Table[i,{i,-20 (π/180),20 (π/180),5 (π/180)}],None}},AspectRatio->1,ImageSize->Round[((72*8)/2.54)],
PlotLegends->BarLegend[{ColorBlend,{-(1.5/100),(1.5/100)}},LegendLabel->Text[Style["τ",FontFamily->"Palatino Linotype",11]],LegendMarkerSize->Round[(72/2.54)7]],LabelStyle->Directive[FontFamily->"Palatino Linotype",Black,11]]

Export["LogicaControlDens.pdf",Rasterize[FiguraD,RasterSize->8000]]

Plot3D[νac,{θ[t],-1.25 θfin,1.25 θfin},{Θ[t],-1.25 ωfin,1.25 ωfin},PlotRange->Full,PlotPoints->200,Mesh->20,MeshStyle->Directive[{Dashed,Borde}],ColorFunction->Table[Blend[{RGBColor["#2a6898"],White,RGBColor["#8fb3d1"]},x],{x,0,1,1/10}],Boxed->True,AxesLabel->{θ,θ',τ},AxesStyle->Directive[FontFamily->"Cambria Math",Black,12],AspectRatio->1,ImageSize->UpTo[500],Exclusions->None,BoxRatios->{1,1,0.75},Ticks->{{-0.1,0.1},{-15,-7.5,7.5,15},{-0.75,-0.5,-0.25,0.25,0.5,0.75}},ViewPoint->{π/4,π/2,π/4},PlotTheme->"Scientific"];

Parámetros de Amortiguamiento

Determinación de coeficientes por ensayo V3

El tiempo se mide directamente en el ensayo, asegurando que las mediciones se realizan en el instante que se indica. La correspondencia con el tiempo real fue verificada mediante comparación con cronómetro.

Carga y adecuación de ensayos

MedSeg=20;
(*Mediciones realizadas*)
Medicion1=Import["C:\\Users\\User\\Free Decay Experiments\\1DL.txt","Table"];
Medicion2=Import["C:\\Users\\User\\Free Decay Experiments\\2DL.txt","Table"];Medicion3=Import["C:\\Users\\User\\Free Decay Experiments\\3DL.txt","Table"];Medicion4=Import["C:\\Users\\User\\Free Decay Experiments\\6DL.txt","Table"];Medicion5=Import["C:\\Users\\User\\Free Decay Experiments\\9DL.txt","Table"];Medicion6=Import["C:\\Users\\User\\Free Decay Experiments\\10DL.txt","Table"];Medicion7=Import["C:\\Users\\User\\Free Decay Experiments\\12DL.txt","Table"];Medicion8=Import["C:\\Users\\User\\Free Decay Experiments\\13DL.txt","Table"];Medicion9=Import["C:\\Users\\User\\Free Decay Experiments\\14DL.txt","Table"];Medicion10=Import["C:\\Users\\User\\Free Decay Experiments\\15DL.txt","Table"];
Medicion11=Import["C:\\Users\\User\\Free Decay Experiments\\16DL.txt","Table"];
Medicion12=Import["C:\\Users\\User\\Free Decay Experiments\\17DL.txt","Table"];
Medicion13=Import["C:\\Users\\User\\Free Decay Experiments\\18DL.txt","Table"];
Medicion14=Import["C:\\Users\\User\\Free Decay Experiments\\19DL.txt","Table"];
Medicion15=Import["C:\\Users\\User\\Free Decay Experiments\\20DL.txt","Table"];

m1=61;m2=43.9;m3=23.7;m4=43.9;me=83.7; (*Pesos en g*)
listC={}; (*Lista de los coeficientes de amortiguamientos viscosos obtenidos para cada ensayo*)
listTr={};(*Lista de los coeficientes de amortiguamientos viscosos obtenidos para cada ensayo*)
Do[

Ensayo=ToExpression["Medicion"<>ToString[j],StandardForm];  (*Selección de medicion para estudio*)
Ensayopos=Table[{Ensayo[[i,1]]-Ensayo[[1,1]],((π Ensayo[[i,2]])/180)},{i,1,Length[Ensayo]}];
Ensayovel=Table[{Ensayo[[i,1]]-Ensayo[[1,1]],((Ensayo[[i+1,2]]-Ensayo[[i,2]])/(Ensayo[[2,1]]-Ensayo[[1,1]]))(π/180)},{i,1,Length[Ensayo]-1}];

(*Obtención de los picos de oscilación*)

PicosPos={}; PicosNeg={};
Do[If[Ensayopos[[i-1,2]]>Ensayopos[[i,2]] ∧ Ensayopos[[i+1,2]]>=Ensayopos[[i,2]] ∧ Ensayopos[[i,2]]<-0.045,AppendTo[PicosNeg,{((i-1)/MedSeg),Ensayopos[[i,2]]}]],{i,2,Length[Ensayo]-1}];
Do[If[Ensayopos[[i-1,2]]<Ensayopos[[i,2]] ∧ Ensayopos[[i+1,2]]<=Ensayopos[[i,2]] ∧ Ensayopos[[i,2]]>0.045,AppendTo[PicosPos,{((i-1)/MedSeg),Ensayopos[[i,2]]}]],{i,2,Length[Ensayo]-1}];

(*Selección de picos para estudio*)

AngMax=16; (*Ángulo máximo para determinación de coeficientes*)
PicosNegLimitados={};
PicosPosLimitados={};
Do[If[Abs[PicosNeg[[i,2]]]<(AngMax(π/180)),AppendTo[PicosNegLimitados,PicosNeg[[i]]]],{i,2,Length[PicosNeg]}]; 
Do[If[Abs[PicosPos[[i,2]]]<(AngMax(π/180)),AppendTo[PicosPosLimitados,PicosPos[[i]]]],{i,2,Length[PicosPos]}];

Show[{ListLinePlot[Ensayopos],ListPlot[PicosPosLimitados,PlotStyle->Red],ListPlot[PicosNegLimitados,PlotStyle->Blue]}];

(*Obtención de la frecuencia natural amortiguada*)

ωdNeg={};(*Frecuencias naturales amortiguadas para ciclos negativos*)
ωdPos={};(*Frecuencias naturales amortiguadas para ciclos positivos*)
Do[AppendTo[ωdNeg,((2π)/((PicosNegLimitados[[i+1,1]]-PicosNegLimitados[[i,1]])))],{i,1,Length[PicosNegLimitados]-1}];(*Frecuencia angular*)
Do[AppendTo[ωdPos,((2π)/((PicosPosLimitados[[i+1,1]]-PicosPosLimitados[[i,1]])))],{i,1,Length[PicosPosLimitados]-1}];(*Frecuencia angular*)
ωdprom=((1/(Length[ωdNeg]))Sum[ωdNeg[[i]],{i,1,Length[ωdNeg]}]+(1/(Length[ωdPos]))Sum[ωdPos[[i]],{i,1,Length[ωdPos]}])/2;(* Print["(Overscript[(Subscript[ω, d]), -]): ",ωdprom//N];*)
(*ωgraph=Show[ListLinePlot[{ωdPos,ωdNeg},AxesLabel->{"|(Subscript[θ, k])|",(Subscript[ω, d])[Hz]},AxesStyle->Directive[FontFamily->"Cambria", 14],PlotStyle->{RGBColor["#9C27B0"],RGBColor["#E91E63"]}],Graphics[{Text[Style["(Overscript[(Subscript[ω, d]), -])",FontFamily->"Cambria",14],{0.75,ωdprom+0.02}],Dashed,Line[{{0,ωdprom},{Length[ωdNeg],ωdprom}}]}]];*)

(*Obtención de coeficiente de amortiguamiento viscoso*)

PicosComb=Sort[Join[Abs[PicosNegLimitados],PicosPosLimitados]]; (*Lista de picos combinados*)
δComb={}; (*Lista de decrementos*)
Do[AppendTo[δComb,{Abs[PicosComb[[i,2]]],(10^4)Abs[(1/((i-1)π))Log[((Abs[PicosComb[[1,2]]]-Abs[PicosComb[[2,2]]])/(Abs[PicosComb[[i-1,2]]]-Abs[PicosComb[[i,2]]]))]]}],{i,6,Length@PicosComb}];(*(10^4) para corrección por unidades*)
δpromComb=(1/(Length@δComb))(∑_(i=1)^(Length@δComb))δComb[[i,2]];(*Valor promedio de decremento para el ensayo*)
CComb=Table[{i,δComb[[i,2]]*Mp*ωdprom},{i,1,Length@δComb}];
CpromComb=δpromComb*Mp*ωdprom ;(*Print["c: ",CpromComb];*)(*Coeficiente de amortiguamiento viscoso*)

(*Determinación de los coeficientes de Coulomb según sentido de oscilación*)

(*(10^6) para corrección por unidades*)
TrNeg={}; (*Torque residual para oscilaciones anti-horarias*)
Do[AppendTo[TrNeg,{Abs[PicosNegLimitados[[i,2]]],(10^6)((ωdprom)^2)*((Abs[PicosNegLimitados[[2,2]]]*Abs[PicosNegLimitados[[i-1,2]]]-Abs[PicosNegLimitados[[1,2]]]*Abs[PicosNegLimitados[[i,2]]])/(Abs[PicosNegLimitados[[2,2]]]+Abs[PicosNegLimitados[[1,2]]]-Abs[PicosNegLimitados[[i-1,2]]]-Abs[PicosNegLimitados[[i,2]]]))}],{i,3,Length@PicosNegLimitados}];
TrPos={};(*Torque residual para oscilaciones horarias*)
Do[AppendTo[TrPos,{Abs[PicosPosLimitados[[i,2]]],(10^6)((ωdprom)^2)*((Abs[PicosPosLimitados[[2,2]]]*Abs[PicosPosLimitados[[i-1,2]]]-Abs[PicosPosLimitados[[1,2]]]*Abs[PicosPosLimitados[[i,2]]])/(Abs[PicosPosLimitados[[2,2]]]+Abs[PicosPosLimitados[[1,2]]]-Abs[PicosPosLimitados[[i-1,2]]]-Abs[PicosPosLimitados[[i,2]]]))}],{i,3,Length@PicosPosLimitados}];
TrCCW=(1/(Length@TrNeg))(∑_(i=1)^(Length@TrNeg))TrNeg[[i,2]];(*Print["(Subscript[T, r-]): ",TrCCW];*)(*Coeficiente de fricción seca anti-horario*)
TrCW=(1/(Length@TrPos))(∑_(i=1)^(Length@TrPos))TrPos[[i,2]];(*Print["(Subscript[T, r+]): ",TrCW];*)(*Coeficiente de fricción seca horario*)
(*Determinación del coeficiente de Coulomb promediado*)

Trp=(((TrCCW+TrCW))/2); (*Print["(Subscript[T, r]): ",Trp];*)
AppendTo[listC,CpromComb];
AppendTo[listTr,Trp];
Clear[Ensayo],
{j,1,15}]
Print["c: ",listC]
Print["(Subscript[T, r]): ",listTr]

Valores promediados

Cpromtot=Mean[listC]
Trptot=Mean[listTr]

Figuras

Donde puede ser de utilidad verificar la determinación de los picos de cada ciclo

Show[ListLinePlot[Ensayopos,PlotStyle->{RGBColor["#1E2A50"],Dashing[{0.005,0.005}]},AspectRatio->1/2,ImageSize->UpTo[600],AxesLabel->{t,"ϕ(t)"},AxesStyle->Directive[FontFamily->"Cambria", 14],PlotRange->{{0,Ensayopos[[Length@Ensayopos,1]]},{-1.1π,1.1π}},Ticks->{{10,20,30,40,50,60},{-π,-π/2,-π/4,0,π/4,π/2,π}}],ListPlot[{PicosPos,PicosNeg},PlotStyle->{{RGBColor["#E91E63"],PointSize[0.008]},{RGBColor["#9C27B0"],PointSize[0.008]}}],Graphics[{Dashed,Line[{{0,AngMax(π/180)},{65,AngMax(π/180)}}],Line[{{0,-AngMax(π/180)},{65,-AngMax(π/180)}}]}],AspectRatio->1/GoldenRatio,ImageSize->UpTo[600]]
(*Frecuencia natural amortiguada*)
Show[ListLinePlot[{ωdPos,ωdNeg},AxesLabel->{"|(Subscript[θ, k])|",(Subscript[ω, d])},AxesStyle->Directive[FontFamily->"Cambria", 14],PlotStyle->{RGBColor["#9C27B0"],RGBColor["#E91E63"]}],Graphics[{Text[Style["(Overscript[(Subscript[ω, d]), -])",FontFamily->"Cambria",14],{0.75,ωdprom+0.02}],Dashed,Line[{{0,ωdprom},{Length[ωdNeg],ωdprom}}]}],AspectRatio->1/GoldenRatio,ImageSize->UpTo[400]]
(*Coeficiente de amortiguamiento viscoso*)
Show[ListLinePlot[CComb,PlotRange->Full,PlotStyle->{RGBColor["#1E2A50"],Dashing[{0.005,0.005}]},AxesLabel->{t,"c"},AxesStyle->Directive[FontFamily->"Cambria", 14],AspectRatio->1/GoldenRatio,ImageSize->UpTo[400]],ListPlot[CComb,PlotStyle->{RGBColor["#E91E63"]}]]
(*Coeficientes de amortiguamiento seco*)
Show[ListLinePlot[{TrNeg,TrPos},PlotStyle->Directive[RGBColor["#1E2A50"],Dashing[{0.005,0.005}]],AxesLabel->{t,"(Subscript[T, r-]),(Subscript[T, r+])"},AxesStyle->Directive[FontFamily->"Cambria", 14],AspectRatio->1/GoldenRatio,ImageSize->UpTo[400]],ListPlot[{TrNeg,TrPos},PlotStyle->{RGBColor["#9C27B0"],RGBColor["#E91E63"]}]]

La distribución de los coeficientes es:

listAmort={};
Do[AppendTo[listAmort,{listC[[i]],listTr[[i]]}],{i,1,Length@listC}]
Show[{ListPlot[listAmort,PlotRange->{{0,200000},{200000,400000}},AxesLabel->{"C","(Subscript[T, r])"},AxesStyle->Directive[FontFamily->"Cambria",12],ImageSize->UpTo[500],AspectRatio->1],Graphics[{Red,PointSize[0.02],Point[{Cpromtot,Trptot}]}]}]

Comparativa simulación-experimento

Ensayo=Medicion12;  (*Selección de medicion para estudio: 4 10 12 17 20 21*)
Ensayopos=Table[{Ensayo[[i,1]]-Ensayo[[1,1]],((π Ensayo[[i,2]])/180)},{i,1,Length[Ensayo]}];
Ensayovel=Table[{Ensayo[[i,1]]-Ensayo[[1,1]],((Ensayo[[i+1,2]]-Ensayo[[i,2]])/(Ensayo[[2,1]]-Ensayo[[1,1]]))(π/180)},{i,1,Length[Ensayo]-1}];
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*)
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*)
Icm=(lc^2) (m1+m2+m3+m4); (*Momento de inercia del centro de masa*)
Icorr=330000; (*Momento de inercia de corrección en base a ensayo de I=τ/α, 663490 según ensayos*)
Hex=0; (*Amplitud de forzamiento horizontal*)
Ωex=8;(*Frecuencia del forzamiento horizontal*)
X[t]=Hex Cos[Ωex t];
Hey=0; (*Amplitud de forzamiento vertical*)
Ωey=8; (*Frecuencia del forzamiento vertical*)
Y[t]=Hey Cos[Ωey t];
tmax=80;
EqMovθ=(Imasas+Iz+Icm+Icorr) θ''[t]+Cpromtot  θ'[t]+Trptot Sign[θ'[t]]+(lc(∑_(i=1)^(n-1))m[[i]]D[X[t],{t,2}])Cos[θ[t]]+(g(l1 m1-l3 m3)+lc(∑_(i=1)^(n-1))m[[i]] D[Y[t],{t,2}])Sin[θ[t]]==0;
SolMovθ=NDSolve[{EqMovθ,θ[0]==Ensayo[[1,2]](π/180),θ'[0]==((Ensayo[[2,2]]-Ensayo[[1,2]])/(Ensayo[[2,1]]-Ensayo[[1,1]]))(π/180)},{θ,θ'},{t,0,tmax},MaxStepSize->Infinity,AccuracyGoal->8];

Ensayovel=Table[{((i-1)/MedSeg),6π(Ensayopos[[i,2]]-Ensayopos[[i-1,2]])},{i,2,Length@Ensayopos-1}];(*Revisar porque sale con 6π*)
EspacioParametros=Table[{Ensayopos[[i,2]],Ensayovel[[i,2]]},{i,1,Length@Ensayovel}];

Para la determinación del espacio de parámetros, es necesario determinar la curva de velocidad del ensayo.

Show[
ListLinePlot[{Ensayopos},PlotRange->{{0,Ensayopos[[Length@Ensayopos,1]]},{-1.1π,1.1π}},PlotStyle->RGBColor["#E91E63"],Ticks->{Automatic,{-π,-π/2,-π/4,π/4,π/2,π}},AspectRatio->1/2,ImageSize->UpTo[800]],
Plot[θ[t]/.SolMovθ,{t,0,tmax},PlotRange->Full,PlotPoints->100,PlotStyle->{Thickness[0.0015],Dashing[{0.002,0.002}],RGBColor["#1E2A50"]}]]

Show[ParametricPlot[{θ[t],θ'[t]}/.SolMovθ,{t,0,tmax},PlotStyle->{RGBColor["#1E2A50"],Thickness[0.003],Dashing[{0.01,0.01}]},PlotRange->{{-π,π},{-3.2π,3.2π}},AxesLabel->{θ[t],θ'[t]},AxesStyle->Directive[FontFamily->"Cambria",14],Ticks->{{-π,π},{-3π,-2π,-π,π,2π,3π}},ImageSize->UpTo[400],AspectRatio->1],
	ListLinePlot[EspacioParametros,PlotStyle->{RGBColor["#E91E63"],Thickness[0.0015]},PlotRange->{{-π,π},Automatic},AxesLabel->{θ[t],θ'[t]},AxesStyle->Directive[FontFamily->"Cambria",14],Ticks->{{-π,π},{-3π,-2π,-π,π,2π,3π}}]]

Clear[Ensayo,Ensayopos,Ensayovel]

Obtención de coeficientes optimizados

A partir de los coeficientes obtenidos mediante el método de Horton, ahora se realiza un barrido en las cercanías de estos parámetros y se buscá el conjunto local (c,(Subscript[T, r])) que minimice el error

cuadrático medio del modelo. Para la determinación del error se utilizan las tablas de Ensayopos y Ensayovel de todos los ensayos realizados.

La rutina desarrollada asigna valores de (c,(Subscript[T, r])) alrededor de los valores obtenidos mediante Horton, y resuelve la ecuación diferencial para luego calcular error cuadrático medio como ((((((Subscript[θ, exp])-(Subscript[θ, sim]))/(Subscript[θ, exp]))))^2). Los valores obtenidos para cada ensayo son sumados y se obtiene un “puntaje” de cercanía. La combinación que menor puntaje obtenga será aquella que mejor haya ajustado los parámetros.

Como se encuentran diferencias en la frecuencia natural amortiguada, se analiza que a la hora de determinar la densidad de la estructura se asumió una densidad constante en todo su volumen, condición que puede no ser verdadera debido a las técnicas de fabricación utilizadas. Sumado a esto y debido a la baja inercia del sistema físico, se evalúa la posibilidad de la modificación del momento de inercia debido a aquel del encoder, eje y rodamientos utilizados en el prototipo fabricado. Debido a estas incertidumbres en cuanto a la medición, se entiende que existen errores en la determinación del momento de inercia debidas a la simplificación del modelo, y a las limitadas herramientas disponibles para verificarlo. En base a esto, el siguiente código busca evaluar el momento de inercia de corrección que más acerque las oscilaciones simuladas a aquellas experimentales.

Para cada ensayo con los coeficientes de Horton ya determinados, se evalúa un rango de estos para distintos valores de momentos de inercia, y se anota el puntaje obtenido. Luego se evaluan los puntajes para todos los conjuntos (c,(Subscript[T, r]),(Subscript[I, corr])), y se elige el de menor puntaje. Se eligen saltos discretos de 10000 (gmm^2) para la determinación de la inercia, y se hace un barrido entre los valores 460000 y 700000 (el mejor ajuste a ojo se da en las cercanías a 500000, pero el ensayo de determinación de la inercia indica que (Subscript[I, corr]) debería ser de 660000). Se evalúan entonces 24 valores de inercia para cada combinación de parámetros de amortiguamiento.

Carga de datos

Nred=2; (*Intervalo cada lo cual se consideran los puntos de las respuestas*)

Lista1={};
Medicion=Medicion1;
Do[AppendTo[Lista1,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista2={};
Medicion=Medicion2;
Do[AppendTo[Lista2,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista3={};
Medicion=Medicion3;
Do[AppendTo[Lista3,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista4={};
Medicion=Medicion4;
Do[AppendTo[Lista4,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista5={};
Medicion=Medicion5;
Do[AppendTo[Lista5,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista6={};
Medicion=Medicion6;
Do[AppendTo[Lista6,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista7={};
Medicion=Medicion7;
Do[AppendTo[Lista7,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista8={};
Medicion=Medicion8; (*Muy mala medición, se elimina de la evaluación*)
Do[AppendTo[Lista8,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista9={};
Medicion=Medicion9;
Do[AppendTo[Lista9,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista10={};
Medicion=Medicion10;
Do[AppendTo[Lista10,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista11={};
Medicion=Medicion11;
Do[AppendTo[Lista11,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista12={};
Medicion=Medicion12;
Do[AppendTo[Lista12,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista13={};
Medicion=Medicion13;
Do[AppendTo[Lista13,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista14={};
Medicion=Medicion14;
Do[AppendTo[Lista14,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]
Lista15={};
Medicion=Medicion15;
Do[AppendTo[Lista15,{Medicion[[i,1]]-Medicion[[1,1]],((π Medicion[[i,2]])/180),((Medicion[[i+1,2]]-Medicion[[i,2]])/(Medicion[[2,1]]-Medicion[[1,1]]))(π/180)}],{i,1,Length@Medicion-1,Nred}];Clear[Medicion]

EnList={Lista1,Lista2,Lista3,Lista4,Lista5,Lista6,Lista7,Lista8,Lista9,Lista10,Lista11,Lista12,Lista13,Lista14,Lista15}; (*Lista con todas las listas de ensayos en estudio*)

Módulo de optimización

g=9810; (*mm/(s^2)*)
l1=85;l2=85;l3=85;l4=85; (*Longitudes en mm*)
m1=61;m2=43.9;m3=23.7;m4=43.9;me=83.7; (*Pesos en g*)
Iz=347125; (*Error en (gmm^2) debido a diferencias de densidad por la impresión 3D, y que no se considera al encoder, eje y rodamientos en el cálculo*)
Imasas=krel[[1]] m1+krel[[2]] m2+krel[[3]] m3+krel[[4]] m4;
Icm=(lc^2) (m1+m2+m3+m4);
Icorr=327500; (*Momento de inercia de corrección en base a ensayo de I=τ/α. 663490 por cálculo*)
X[t]=0; (*Anulado para ensayo libre*)
Hex=0;
Ωex=0;
Y[t]=0;(*Anulado para ensayo libre*)
Clear[Icorr]
FunErr[Trsim_,Csim_,Icorr_]:=Module[{t0,t1,θ0,v0,SolMovθ,θsolu},
TabAj={};
Do[
(*Resolución de la ecuación diferencial*)
t0=EnList[[k,1,1]]; (*Tiempo inicial*)
t1=EnList[[k,Length@EnList[[k]],1]]; (*Tiempo final*)
θ0=EnList[[k,1,2]];(*Posición inicial*)
v0=EnList[[k,1,3]]; (*Velocidad inicial*)
SolMovθ=NDSolve[{Θ[t]==θ'[t],(Imasas+Iz+Icm+Icorr)Θ'[t]+Csim Θ[t]+Trsim Sign[Θ[t]]+(lc(∑_(i=1)^(n-1))m[[i]]D[X[t],{t,2}])Cos[θ[t]]+(g(l1 m1-l3 m3)+lc(∑_(i=1)^(n-1))m[[i]] D[Y[t],{t,2}])Sin[θ[t]]==0,θ[0]==θ0,Θ[0]==v0},{θ[t],Θ[t]},{t,t0,t1}];
θsolu=θ[t]/.SolMovθ;
Θsolu=Θ[t]/.SolMovθ;
Do[{
(*Evaluación de la solución para cada instancia temporal del ensayo*)
Which[Abs[EnList[[k,i,2]]]>(π/18),AppendTo[TabAj,(((((EnList[[k,i,2]]-(θsolu/.t->EnList[[k,i,1]]))/(EnList[[k,i,2]]))))^2)[[1]]]]},{i,1,Length@EnList[[k]]}]
,{k,1,Length@EnList}];
Total[TabAj]]

FunErr[Trptot,Cpromtot,100000]

Optimización

Evaluación en paralelo de todas las combinaciones ((Subscript[T, r]),c,(Subscript[I, corr])) y generación de tabla de resultados ((Subscript[T, r]),c,(Subscript[I, corr]),Error). Se ajusta en primer lugar en función del momento de inercia Icorr, y luego se deben ajustar los coeficientes de amortiguamiento.

$DefaultKernels[6];
AbsoluteTiming[RegErr={};
SetSharedVariable[RegErr];
ErrCuadRegion=Quiet[ParallelDo[AppendTo[RegErr,{i,j,k,FunErr[i,j,k]}],{i,Trptot,Trptot},{j,Cpromtot,Cpromtot},{k,0,1000000,10000}]]]

ListaIcorr={};
Do[AppendTo[ListaIcorr,{RegErr[[i,3]],RegErr[[i,4]]}],{i,1,Length@RegErr}]
ListLinePlot[ListaIcorr]
ListaIcorr=SortBy[ListaIcorr,Last]; ListaIcorr[[1]]

Se realiza ahora la optimización de los coeficientes de amortiguamiento en una región cercana al punto mínimo del error para el momento de inercia.

Trptot*(10^(-9))
Cpromtot*(10^(-9))

$DefaultKernels[6];
AbsoluteTiming[RegErr={};
SetSharedVariable[RegErr];
ErrCuadRegion=ParallelDo[AppendTo[RegErr,{i*(10^(-9)),j*(10^(-9)),k,FunErr[i,j,k]}],{i,Trptot-80000,Trptot+120000,(200000/30)},{j,Cpromtot-80000,Cpromtot+120000,(200000/30)},{k,700000,900000,15000},Method->"FinestGrained"]]

Guardado de datos

Export["CoeficientesOptimizadosFINAL.txt",RegErr,"Table"]

Carga de datos

RegErr=Import["D:\\Onedrive\\Beca Doctoral CONICET\\Papers y Congresos\\2025 - Rotation control of a parametrically forced magnetic pendulum\\Ensayos\\CoeficientesOptimizados.txt","Table"];

RegErr=Import["C:\\Users\\Nicolas\\OneDrive\\Beca Doctoral CONICET\\Papers y Congresos\\2025 - Rotation control of a parametrically forced magnetic pendulum\\Ensayos\\CoeficientesOptimizados.txt","Table"];

Obtención del conjunto de menor error

MinimalBy[RegErr,Last](*Busca en la lista completa el menor valor del error*)
MaximalBy[RegErr,Last](*Busca en la lista completa el menor valor del error*)

Ordenamiento de listas según momento de inercia

ListasI=GatherBy[RegErr,#[[3]]&]; (*Agrupa los elementos de RegErr en sublistas para inercias de corrección iguales*)

Gráfico de los puntos calculados para cada valor de inercia

ListasI[[5,1]]

ListasIReducidas1 = Map[#[[1 ;; 2]]~Append~#[[4]]&,ListasI[[5]],{1} ];
Minimo440=MinimalBy[ListasIReducidas1 ,Last];
Position[ListasIReducidas1,First[Minimo440]];
Paletaz[z_]:=RGBColor[0.5z,0.65-z,1]; (*Paleta de colores para superficie*)
Show[ListPlot3D[{ListasIReducidas1},PlotRange->{{230000,350000},{55000,125000}},ColorFunction->(Paletaz[#3]&),PlotStyle->Directive[PointSize[0.0175],Opacity[0.5]],ImageSize->Round[72*(12/2.54)],PlotLegends->BarLegend["Paletaz",{0,40000}],Boxed->False,AxesOrigin->Automatic,ViewPoint->{2,0,-0.1 }],
Graphics3D[{PointSize[0.02],Opacity[1],RGBColor["#303F9F"],Point[{Minimo440[[1]]}],RGBColor["#E91E63"],Point[{Trptot,Cpromtot,FunErr[Trptot,Cpromtot,ListasI[[1,1,3]]]}]}]];

ListasI[[5,1,3]]
Minimo440

Gráfico como superficie de contorno

ListaPuntos={};
Do[AppendTo[ListaPuntos,{((ListasIReducidas1[[i,1]])/(10^9)),((ListasIReducidas1[[i,2]])/(10^9)),ListasIReducidas1[[i,3]]}],{i,1,Length@ListasIReducidas1}]

ColorBlend[z_]:=Blend[{{0,RGBColor["#527fba"]},{20000,RGBColor["#E1F5FE"]}},z];
DampingSpaceContour=Show[{ListContourPlot[ListaPuntos,PlotRange->{{0.00025,0.0004},{0.000045,0.00015}},ColorFunction->(ColorBlend[#1]&),ColorFunctionScaling->False,FrameTicks->{{Automatic,None},{Automatic,None}},Contours->{3500,5000,7500,10000,12500,17000,20000},ContourStyle->Directive[Dashed,Opacity[0.25]],AspectRatio->1,ImageSize->Round[72*(8/2.54)],PlotLegends->Placed[Automatic,Below],LabelStyle->Directive[FontFamily->"Palatino Linotype",Black,12],FrameStyle->Directive[FontFamily->"Palatino Linotype", 12]],
Graphics[{Text[Style["(Subscript[P, Horton])",White,FontFamily->"Palatino Linotype",12],{Trptot,105000}],Text[Style["(Subscript[P, MMSE])",White,FontFamily->"Palatino Linotype",12],{Minimo440[[1,1]],72000}],PointSize[0.025],Opacity[1],RGBColor["#1A237E"],Point[{Minimo440[[1,1]],Minimo440[[1,2]]}],RGBColor["#B71C1C"],Point[{Trptot,Cpromtot}]}]}]

ListasIReducidas1[[1]]

Export["DampingSpaceContour.pdf",DampingSpaceContour]

tamaño=8;
DampingSpaceDensity=Show[{ListDensityPlot[ListaPuntos,PlotRange->{{0.000235,0.00035},{0.00006,0.000125}},ColorFunction->(ColorBlend[#1]&),ColorFunctionScaling->False,AspectRatio->1,ImageSize->Round[72*((tamaño)/2.54)],FrameLabel->{"(Subscript[T, r])","c"},FrameTicks->{{Table[i,{i,0,0.0002,0.00002}],None},{Table[i,{i,0,0.0005,0.00005}],None}},FrameStyle->Directive[FontFamily->"Palatino Linotype", 11,Thickness[0.005]]],
Graphics[{Text[Style["(Subscript[P, H])",Bold,Black,FontFamily->"Palatino Linotype",11],{(Trptot/(10^9)),(104000/(10^9))}],Text[Style["(Subscript[P, M])",Bold,Black,FontFamily->"Palatino Linotype",11],{((Minimo440[[1,1]])/(10^9)),(71000/(10^9))}],PointSize[0.03],Opacity[1],RGBColor["#1A237E"],Point[{((Minimo440[[1,1]])/(10^9)),((Minimo440[[1,2]])/(10^9))}],RGBColor["#B71C1C"],Point[{(Trptot/(10^9)),(Cpromtot/(10^9))}]}]}];
DampingBar=BarLegend[{ColorBlend,{0,20000}},LegendLabel->Text[Style["(Subscript[S, MMSE])((Subscript[T, r]),c)",FontFamily->"Palatino Linotype",10]],LabelStyle->Directive[FontFamily->"Palatino Linotype",11],LegendMarkerSize->Round[(72/2.54)tamaño*0.7]];
DampingSpaceDensity=GraphicsRow[{DampingSpaceDensity,DampingBar}]

Export["DampingSpaceDensity.pdf",Rasterize[DampingSpaceDensity,RasterSize->8000]]

Comparativa Horton-Iterativo

Gráfico Horton

g=9810; (*mm/(s^2)*)
l1=85;l2=45;l3=65;l4=45; (*Longitudes en mm*)
m1=61;m2=13.3;m3=23.7;m4=13.3;me=66.7; (*Pesos en g*)
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*)
Icm=(lc^2) (m1+m2+m3+m4); (*Momento de inercia del centro de masa*)
Hey=0; (*Amplitud de forzamiento vertical*)
Ωey=8; (*Frecuencia del forzamiento vertical*)
Y[t]=Hey Cos[Ωey t];
tmax=80;
a=10; (*Ensayo a comparar*)
θ0=EnList[[a,1,2]];
v0=((EnList[[a,2,2]]-EnList[[a,1,2]])/(EnList[[a,2,1]]-EnList[[a,1,1]]));
SolMovθ=NDSolve[{
θ'[t]==Θ[t],
(Imasas+Iz+Icm+ListasI[[a,1,3]]*3) Θ'[t]+Cpromtot Θ[t]+Trptot Sign[Θ[t]]+(lc(∑_(i=1)^(n-1))m[[i]]D[X[t],{t,2}])Cos[θ[t]]+(g(l1 m1-l3 m3)+lc(∑_(i=1)^(n-1))m[[i]] D[Y[t],{t,2}])Sin[θ[t]]==0,θ[0]==θ0,Θ[0]==v0},{θ,Θ},{t,0,EnList[[a,Length@EnList[[a]],1]]}] ; θsolu=θ[t]/.SolMovθ;

GraphHorton=Show[
Plot[θsolu,{t,0,EnList[[a,Length@EnList[[a]],1]]},PlotPoints->300,PlotRange->{{0,EnList[[a,Length@EnList[[a]],1]]},{-π,π}},
PlotStyle->{Thickness[0.0045],RGBColor["#3949AB"]},LabelStyle->{FontFamily->"Palatino Linotype",Black,11},
Axes->False,Frame->True,FrameLabel->{None,"θ(t)"},FrameTicks->{{Table[i,{i,-π,π,π/2}],None},{Table[i,{i,0,50,10}],None}},FrameStyle->Directive[Thickness[0.004]],PlotLabels->{Callout["Simulation",{35,-(π/4)}]},GridLinesStyle->Directive[RGBColor["#757575"], Dashed],GridLines->{ {0},{0}},AspectRatio->(3/4),ImageSize->Round[((72*9)/2.54)]],

ListPlot[EnList[[a]][[All,1;;2]],PlotRange->Full,Joined->True,PlotStyle->{RGBColor["#D81B60"],Thickness[0.0025],Dashing[{0.005,0.005}]},PlotMarkers->{"OpenMarkers",2},PlotLabels->{Callout["Experimental",{35,(π/4)},LabelStyle->{FontFamily->"Palatino Linotype", 11,RGBColor["#D81B60"]}]}],Epilog->{Text[Style["a",FontFamily->"Palatino Linotype",11],{50,-1.6π/2}]}];

Gráfico iterativo

SolMovθ=NDSolve[{
θ'[t]==Θ[t],
(Imasas+Iz+Icm+ListasI[[a,1,3]]*3) Θ'[t]+Minimo440[[1,2]]Θ[t]+Minimo440[[1,1]] Sign[Θ[t]]+(g(l1 m1-l3 m3)+lc(∑_(i=1)^(n-1))m[[i]] D[Y[t],{t,2}])Sin[θ[t]]==0,
θ[0]==θ0,Θ[0]==v0},{θ,Θ},{t,0,EnList[[a,Length@EnList[[a]],1]]}];
θsolu=θ[t]/.SolMovθ;

GraphIterativo=Show[
Plot[θsolu,{t,0,EnList[[a,Length@EnList[[a]],1]]},PlotPoints->300,PlotRange->{{0,EnList[[a,Length@EnList[[a]],1]]},{-π,π}},
PlotStyle->{Thickness[0.0045],RGBColor["#3949AB"]},LabelStyle->{FontFamily->"Palatino Linotype",Black,11},
Axes->False,Frame->True,FrameLabel->{None,"θ(t)"},FrameTicks->{{Table[i,{i,-π,π,π/2}],None},{Table[i,{i,0,50,10}],None}},FrameStyle->Directive[Thickness[0.004]],PlotLabels->{Callout["Simulation",{35,-(π/4)}]},GridLinesStyle->Directive[RGBColor["#757575"], Dashed],GridLines->{ {0},{0}},AspectRatio->(3/4),ImageSize->Round[((72*9)/2.54)]],

ListPlot[EnList[[a]][[All,1;;2]],PlotRange->Full,Joined->True,PlotStyle->{RGBColor["#D81B60"],Thickness[0.0025],Dashing[{0.005,0.005}]},PlotMarkers->{"OpenMarkers",2},PlotLabels->{Callout["Experimental",{35,(π/4)},LabelStyle->{FontFamily->"Palatino Linotype", 11,RGBColor["#D81B60"]}]}],Epilog->{Text[Style["b",FontFamily->"Palatino Linotype",11],{50,-1.6π/2}]}];

6

Combinación de gráficos

FiguraDecr=GraphicsColumn[{GraphHorton,GraphIterativo}]

Export["Horton-Iterative.pdf",FiguraDecr]

Verificación de ajuste. Repetir gráfico de espacio de parámetros para los coeficientes ajustados comparados con los datos experimentales.

Ensayopos=Table[{EnList[[a,i,1]],EnList[[a,i,2]]},{i,1,Length[EnList[[a]]]}];
Ensayovel=Table[{EnList[[a,i,1]],EnList[[a,i,3]]},{i,1,Length[EnList[[a]]]}];
EspacioParametros=Table[{Ensayopos[[i,2]],Ensayovel[[i,2]]},{i,1,Length[EnList[[a]]]}];

Show[
ListLinePlot[{Ensayopos},PlotRange->{{0,Ensayopos[[Length@Ensayopos,1]]},{-1.1π,1.1π}},PlotStyle->RGBColor["#E91E63"],Ticks->{Automatic,{-π,-π/2,-π/4,π/4,π/2,π}},AspectRatio->1/2,ImageSize->UpTo[800]],
Plot[θ[t]/.SolMovθ,{t,0,tmax},PlotRange->Full,PlotPoints->100,PlotStyle->{Thickness[0.0015],Dashing[{0.002,0.002}],RGBColor["#1E2A50"]}]]

Show[ParametricPlot[{θ[t],θ'[t]}/.SolMovθ,{t,0,tmax},PlotStyle->{RGBColor["#1E2A50"],Thickness[0.003],Dashing[{0.01,0.01}]},PlotRange->{{-π,π},{-3.2π,3.2π}},AxesLabel->{θ[t],θ'[t]},AxesStyle->Directive[FontFamily->"Cambria",14],Ticks->{{-π,π},{-3π,-2π,-π,π,2π,3π}},ImageSize->UpTo[400],AspectRatio->1],
	ListLinePlot[EspacioParametros,PlotStyle->{RGBColor["#E91E63"],Thickness[0.0015]},PlotRange->{{-π,π},Full},AxesLabel->{θ[t],θ'[t]},AxesStyle->Directive[FontFamily->"Cambria",14],Ticks->{{-π,π},{-3π,-2π,-π,π,2π,3π}}]]