AJO_1415_Ejemplo_poligonal_Tindouf

940 days ago by sevillad1415

Dados tres puntos A, B, C conocemos las coordenadas de A, las tres distancias y los tres acimutes, con las precisiones indicadas.

Debajo está el ajuste con las ocho observaciones (caso 1), descartando una observación (casos 2a y 2b), y descartando dos observaciones (caso 3). La hoja está preparada para ser ejecutada de principio a fin (menú "Action" / "Evaluate all"). Se comparan los cuatro métodos mostrando juntas las elipses de error.

 

Caso 1: con las ocho mediciones

reset() 
       
lxA,lyA,sxA,syA = 711708,2969077,1,1 ldAB,sdAB = 120513.825,0.050 laAB,saAB = -58.8554,0.0020 ldBC,sdBC = 33309.819,0.029 laBC,saBC = -38.5767,0.0020 ldAC,sdAC = 152504.929,0.057 laAC,saAC = -54.4977,0.0020 SigmaL = diagonal_matrix([sxA^2,syA^2,sdAB^2,saAB^2,sdBC^2,saBC^2,sdAC^2,saAC^2]) QL = SigmaL P = QL^(-1) var('xA,yA,xB,yB,xC,yC') F(xA,yA,xB,yB,xC,yC) = ( xA , yA , sqrt((xA-xB)^2+(yA-yB)^2) , 200/pi*atan2(xB-xA,yB-yA) , sqrt((xB-xC)^2+(yB-yC)^2) , 200/pi*atan2(xC-xB,yC-yB) , sqrt((xA-xC)^2+(yA-yC)^2) , 200/pi*atan2(xC-xA,yC-yA) ) 
       
lxA+ldAB*sin(pi/200*laAB).n(),lyA+ldAB*cos(pi/200*laAB).n(); lxA+ldAC*sin(pi/200*laAC).n(),lyA+ldAC*cos(pi/200*laAC).n() 
       


X0 = vector([711708,2969077,615499,3.04165*10^6,596527,3.069032*10^6]).column() F(*X0.column(0)).n() 
       
point((X0[0,0],X0[1,0]),figsize=4,aspect_ratio=1) + point((X0[2,0],X0[3,0])) + point((X0[4,0],X0[5,0])) 
       
L = vector([lxA,lyA,ldAB,laAB,ldBC,laBC,ldAC,laAC]).column() L 
       
J0 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X0.column(0)).n() N0 = transpose(J0)*P*J0 QX0 = N0^(-1) SigmaX0 = QX0 SigmaX0; map(sqrt,SigmaX0.diagonal()) 
       


V0 = L - F(*X0.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] sumaV0 
       
DeltaL0 = L - F(*X0.column(0)).column() DeltaX0 = (transpose(J0)*P*J0)^(-1) * (transpose(J0)*P*DeltaL0) DeltaX0 = DeltaX0.n() DeltaX0 
       
X1 = X0 + DeltaX0 V1 = L - F(*X1.column(0)).column().n() sumaV1 = transpose(V1)*P*V1 sumaV1 = sumaV1[0,0] sumaV0; sumaV1 
       


J1 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X1.column(0)).n() DeltaL1 = L - F(*X1.column(0)).column() DeltaX1 = (transpose(J1)*P*J1)^(-1) * (transpose(J1)*P*DeltaL1) DeltaX1 = DeltaX1.n() DeltaX1 
       
X2 = X1 + DeltaX1 V2 = L - F(*X2.column(0)).column().n() sumaV2 = transpose(V2)*P*V2 sumaV2 = sumaV2[0,0] sumaV0; sumaV1; sumaV2 
       




X2 
       
J2 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X2.column(0)).n() N2 = transpose(J2)*P*J2 QX2 = N2^(-1) SigmaX2 = QX2 SigmaX2; map(sqrt,SigmaX2.diagonal()) 
       


map(sqrt,SigmaX0.diagonal()) 
       
redund = L.nrows() - X0.nrows() vdr_estim = sumaV2/redund vdr_estim 
       
puntoA = X2[0,0],X2[1,0] puntoB = X2[2,0],X2[3,0] puntoC = X2[4,0],X2[5,0] SigmaA = SigmaX2.submatrix(row=0, col=0, nrows=2, ncols=2) SigmaB = SigmaX2.submatrix(row=2, col=2, nrows=2, ncols=2) SigmaC = SigmaX2.submatrix(row=4, col=4, nrows=2, ncols=2) color = "blue" c = 2.447 
       
e = SigmaA.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaA[0,1]/(max(e)-SigmaA[0,0])) plotAcaso1 = ellipse(puntoA, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoA,color=color) plotAcaso1.show(figsize=2) 
       
/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/all_notebook.\
py:3: UserWarning: Using generic algorithm for an inexact ring, which
will probably give incorrect results due to numerical precision issues.
  
/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/all_notebook.py:3: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.
  
e = SigmaB.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaB[0,1]/(max(e)-SigmaB[0,0])) plotBcaso1 = ellipse(puntoB, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoB,color=color) plotBcaso1.show(figsize=2) 
       
e = SigmaC.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaC[0,1]/(max(e)-SigmaC[0,0])) plotCcaso1 = ellipse(puntoC, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoC,color=color) plotCcaso1.show(figsize=2) 
       
 
       

Caso 2a: con siete mediciones (quitamos el acimut de AC)

lxA,lyA,sxA,syA = 711708,2969077,1,1 ldAB,sdAB = 120513.825,0.050 laAB,saAB = -58.8554,0.0020 ldBC,sdBC = 33309.819,0.029 laBC,saBC = -38.5767,0.0020 ldAC,sdAC = 152504.929,0.057 SigmaL = diagonal_matrix([sxA^2,syA^2,sdAB^2,saAB^2,sdBC^2,saBC^2,sdAC^2]) QL = SigmaL P = QL^(-1) var('xA,yA,xB,yB,xC,yC') F(xA,yA,xB,yB,xC,yC) = ( xA , yA , sqrt((xA-xB)^2+(yA-yB)^2) , 200/pi*atan2(xB-xA,yB-yA) , sqrt((xB-xC)^2+(yB-yC)^2) , 200/pi*atan2(xC-xB,yC-yB) , sqrt((xA-xC)^2+(yA-yC)^2) ) L = vector([lxA,lyA,ldAB,laAB,ldBC,laBC,ldAC]).column() L 
       
F(*X0.column(0)).n() 
       
J0 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X0.column(0)).n() N0 = transpose(J0)*P*J0 QX0 = N0^(-1) SigmaX0 = QX0 SigmaX0; map(sqrt,SigmaX0.diagonal()) 
       


V0 = L - F(*X0.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] sumaV0 
       
DeltaL0 = L - F(*X0.column(0)).column() DeltaX0 = (transpose(J0)*P*J0)^(-1) * (transpose(J0)*P*DeltaL0) DeltaX0 = DeltaX0.n() X1 = X0 + DeltaX0 V1 = L - F(*X1.column(0)).column().n() sumaV1 = transpose(V1)*P*V1 sumaV1 = sumaV1[0,0] DeltaX0; sumaV0; sumaV1 
       




J1 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X1.column(0)).n() DeltaL1 = L - F(*X1.column(0)).column() DeltaX1 = (transpose(J1)*P*J1)^(-1) * (transpose(J1)*P*DeltaL1) DeltaX1 = DeltaX1.n() X2 = X1 + DeltaX1 V2 = L - F(*X2.column(0)).column().n() sumaV2 = transpose(V2)*P*V2 sumaV2 = sumaV2[0,0] DeltaX1; sumaV0; sumaV1; sumaV2 
       






X2 
       
J2 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X2.column(0)).n() N2 = transpose(J2)*P*J2 QX2 = N2^(-1) SigmaX2 = QX2 SigmaX2; map(sqrt,SigmaX2.diagonal()) 
       


map(sqrt,SigmaX0.diagonal()) 
       
redund = L.nrows() - X0.nrows() vdr_estim = sumaV2/redund vdr_estim 
       
puntoA = X2[0,0],X2[1,0] puntoB = X2[2,0],X2[3,0] puntoC = X2[4,0],X2[5,0] SigmaA = SigmaX2.submatrix(row=0, col=0, nrows=2, ncols=2) SigmaB = SigmaX2.submatrix(row=2, col=2, nrows=2, ncols=2) SigmaC = SigmaX2.submatrix(row=4, col=4, nrows=2, ncols=2) color = "green" c = 2.447 
       
e = SigmaA.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaA[0,1]/(max(e)-SigmaA[0,0])) plotAcaso2a = ellipse(puntoA, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoA,color=color) plotAcaso2a.show(figsize=2) 
       
e = SigmaB.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaB[0,1]/(max(e)-SigmaB[0,0])) plotBcaso2a = ellipse(puntoB, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoB,color=color) plotBcaso2a.show(figsize=2) 
       
e = SigmaC.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaC[0,1]/(max(e)-SigmaC[0,0])) plotCcaso2a = ellipse(puntoC, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoC,color=color) plotCcaso2a.show(figsize=2) 
       
(plotAcaso1 + plotAcaso2a).show(figsize=4);(plotBcaso1 + plotBcaso2a).show(figsize=4);(plotCcaso1 + plotCcaso2a).show(figsize=4) 
       








Caso 2b: con siete mediciones (quitamos la distancia AC)

lxA,lyA,sxA,syA = 711708,2969077,1,1 ldAB,sdAB = 120513.825,0.050 laAB,saAB = -58.8554,0.0020 ldBC,sdBC = 33309.819,0.029 laBC,saBC = -38.5767,0.0020 laAC,saAC = -54.4977,0.0020 SigmaL = diagonal_matrix([sxA^2,syA^2,sdAB^2,saAB^2,sdBC^2,saBC^2,saAC^2]) QL = SigmaL P = QL^(-1) var('xA,yA,xB,yB,xC,yC') F(xA,yA,xB,yB,xC,yC) = ( xA , yA , sqrt((xA-xB)^2+(yA-yB)^2) , 200/pi*atan2(xB-xA,yB-yA) , sqrt((xB-xC)^2+(yB-yC)^2) , 200/pi*atan2(xC-xB,yC-yB) , 200/pi*atan2(xC-xA,yC-yA) ) L = vector([lxA,lyA,ldAB,laAB,ldBC,laBC,laAC]).column() L 
       
F(*X0.column(0)).n() 
       
J0 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X0.column(0)).n() N0 = transpose(J0)*P*J0 QX0 = N0^(-1) SigmaX0 = QX0 SigmaX0; map(sqrt,SigmaX0.diagonal()) 
       


V0 = L - F(*X0.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] sumaV0 
       
DeltaL0 = L - F(*X0.column(0)).column() DeltaX0 = (transpose(J0)*P*J0)^(-1) * (transpose(J0)*P*DeltaL0) DeltaX0 = DeltaX0.n() X1 = X0 + DeltaX0 V1 = L - F(*X1.column(0)).column().n() sumaV1 = transpose(V1)*P*V1 sumaV1 = sumaV1[0,0] DeltaX0; sumaV0; sumaV1 
       




J1 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X1.column(0)).n() DeltaL1 = L - F(*X1.column(0)).column() DeltaX1 = (transpose(J1)*P*J1)^(-1) * (transpose(J1)*P*DeltaL1) DeltaX1 = DeltaX1.n() X2 = X1 + DeltaX1 V2 = L - F(*X2.column(0)).column().n() sumaV2 = transpose(V2)*P*V2 sumaV2 = sumaV2[0,0] DeltaX1; sumaV0; sumaV1; sumaV2 
       






X2 
       
J2 = jacobian(F,(xA,yA,xB,yB,xC,yC))(*X2.column(0)).n() N2 = transpose(J2)*P*J2 QX2 = N2^(-1) SigmaX2 = QX2 SigmaX2; map(sqrt,SigmaX2.diagonal()) 
       


map(sqrt,SigmaX0.diagonal()) 
       
redund = L.nrows() - X0.nrows() vdr_estim = sumaV2/redund vdr_estim 
       
puntoA = X2[0,0],X2[1,0] puntoB = X2[2,0],X2[3,0] puntoC = X2[4,0],X2[5,0] SigmaA = SigmaX2.submatrix(row=0, col=0, nrows=2, ncols=2) SigmaB = SigmaX2.submatrix(row=2, col=2, nrows=2, ncols=2) SigmaC = SigmaX2.submatrix(row=4, col=4, nrows=2, ncols=2) color = "black" c = 2.447 
       
e = SigmaA.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaA[0,1]/(max(e)-SigmaA[0,0])) plotAcaso2b = ellipse(puntoA, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoA,color=color) plotAcaso2b.show(figsize=2) 
       
e = SigmaB.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaB[0,1]/(max(e)-SigmaB[0,0])) plotBcaso2b = ellipse(puntoB, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoB,color=color) plotBcaso2b.show(figsize=2) 
       
e = SigmaC.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaC[0,1]/(max(e)-SigmaC[0,0])) plotCcaso2b = ellipse(puntoC, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoC,color=color) plotCcaso2b.show(figsize=2) 
       
(plotAcaso1 + plotAcaso2a + plotAcaso2b).show(figsize=4);(plotBcaso1 + plotBcaso2a + plotBcaso2b).show(figsize=4);(plotCcaso1 + plotCcaso2a + plotCcaso2b).show(figsize=4) 
       








 
       

Caso 3: con seis mediciones (quitamos las dos de AC), ¡ya no es un ajuste!

lxA,lyA,sxA,syA = 711708,2969077,1,1 ldAB,sdAB = 120513.825,0.050 laAB,saAB = -58.8554,0.0020 ldBC,sdBC = 33309.819,0.029 laBC,saBC = -38.5767,0.0020 var('xA,yA,dAB,aAB,dBC,aBC') F(xA,yA,dAB,aAB,dBC,aBC) = ( xA , yA , xA+dAB*sin(pi/200*aAB) , yA+dAB*cos(pi/200*aAB) , xA+dAB*sin(pi/200*aAB)+dBC*sin(pi/200*aBC) , yA+dAB*cos(pi/200*aAB)+dBC*cos(pi/200*aBC) ) 
       
X = F(lxA,lyA,ldAB,laAB,ldBC,laBC).column().n() X 
       
J = jacobian(F,(xA,yA,dAB,aAB,dBC,aBC))(lxA,lyA,ldAB,laAB,ldBC,laBC).n() SigmaL = diagonal_matrix([sxA^2,syA^2,sdAB^2,saAB^2,sdBC^2,saBC^2]) SigmaX = J * SigmaL * transpose(J) SigmaX; map(sqrt,SigmaX.diagonal()) 
       


puntoA = X[0,0],X[1,0] puntoB = X[2,0],X[3,0] puntoC = X[4,0],X[5,0] SigmaA = SigmaX.submatrix(row=0, col=0, nrows=2, ncols=2) SigmaB = SigmaX.submatrix(row=2, col=2, nrows=2, ncols=2) SigmaC = SigmaX.submatrix(row=4, col=4, nrows=2, ncols=2) color = "red" c = 2.447 
       
e = SigmaA.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaA[0,1]/(max(e)-SigmaA[0,0])); alfa=0 plotAcaso3 = ellipse(puntoA, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoA,color=color) plotAcaso3.show(figsize=2) 
       
e = SigmaB.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaB[0,1]/(max(e)-SigmaB[0,0])) plotBcaso3 = ellipse(puntoB, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoB,color=color) plotBcaso3.show(figsize=2) 
       
e = SigmaC.eigenvalues() smax = max(map(sqrt,e)) smin = min(map(sqrt,e)) alfa = arctan(SigmaC[0,1]/(max(e)-SigmaC[0,0])) plotCcaso3 = ellipse(puntoC, c*smax, c*smin, pi/2-alfa, figsize=5, aspect_ratio=1,color=color) + point(puntoC,color=color) plotCcaso3.show(figsize=2) 
       
(plotAcaso1 + plotAcaso2a + plotAcaso2b + plotAcaso3).show(figsize=4);(plotBcaso1 + plotBcaso2a + plotBcaso2b + plotBcaso3).show(figsize=4);(plotCcaso1 + plotCcaso2a + plotCcaso2b + plotCcaso3).show(figsize=4)