Ajuste de observaciones por mínimos cuadrados y elipse de error en MATLAB

Resolución de Ajuste de Observaciones mediante Mínimos Cuadrados

A continuación se presenta el script desarrollado en MATLAB para la resolución de un ejercicio de ajuste de observaciones, centrado en el cálculo de coordenadas y la determinación de elipses de error.

1. Configuración y Datos Iniciales

En esta sección se definen los puntos de control, las distancias observadas y sus respectivas precisiones (sigmas).

% RESOLUCIÓN EJERCICIO 1 - AJUSTE DE OBSERVACIONES (MÍNIMOS CUADRADOS)
clear; clc;

% 1. Datos iniciales
% Coordenadas de los puntos conocidos
P_fijos = [30, 150; 
           10, 120; 
           50, 50];

% Distancias observadas (L) y sus precisiones (sigmas)
D_obs = [125.0; 133.5; 98.6];
sigmas = [0.5; 0.2; 0.2];

% Coordenadas aproximadas del punto P
X0 = [140; 90]; 

% Matriz de pesos (P)
% Se asume varianza de referencia a priori = 1
W = diag(1 ./ (sigmas.^2));

2. Proceso Iterativo de Cálculo

Se implementa un bucle para refinar las coordenadas mediante el Jacobiano de las derivadas parciales y las ecuaciones normales.

% 2. Proceso iterativo
iteraciones = 5;
X = X0;
fprintf('--- PROCESO DE ITERACIÓN ---\n');

for i = 1:iteraciones
    % Distancias calculadas desde las coordenadas actuales a los puntos fijos
    dX = X(1) - P_fijos(:,1);
    dY = X(2) - P_fijos(:,2);
    D_cal = sqrt(dX.^2 + dY.^2);

    % Matriz de diseño A (Jacobiano de las derivadas parciales)
    A = [dX./D_cal, dY./D_cal]; 

    % Vector de términos independientes (Observado - Calculado)
    L = D_obs - D_cal;

    % Ecuaciones normales
    N = A' * W * A;
    U = A' * W * L;

    % Vector de correcciones
    delta = N \ U;

    % Actualización de las coordenadas
    X = X + delta; 

    fprintf('ITERACIÓN %d\n', i);
    fprintf('CORRECCIONES: %.2f metros (dx), %.2f metros (dy)\n', delta(1), delta(2));
    fprintf('x0 = %.4f\ny0 = %.4f\n\n', X(1), X(2));

    % Condición de parada
    if norm(delta) < 1e-5
        break;
    end
end

3. Resultados del Ajuste y Precisión

Evaluación de los residuos, la varianza a posteriori y la matriz de covarianza para determinar la precisión de las coordenadas finales.

% Vector de residuos v = A*delta - L
V = D_cal - D_obs; 
fprintf('--- RESULTADOS DEL AJUSTE ---\n');
fprintf('Residuos v:\n'); disp(V);

n = length(D_obs);
u = length(X);
gl = n - u;

% Varianza de referencia a posteriori
var_post = (V' * W * V) / gl;
fprintf('Varianza de referencia a posteriori: %.4f\n\n', var_post);

Sx = inv(N);
fprintf('Matriz Covarianza Sx:\n'); disp(Sx);

% Precisiones de las coordenadas
Sx_x = sqrt(Sx(1,1));
Sx_y = sqrt(Sx(2,2));
fprintf('PRECISIONES DE LAS COORDENADAS SX=%.2f m SY=%.2f m\n\n', Sx_x, Sx_y);

4. Cálculo de la Elipse de Error Estándar

Determinación de los semiejes y la orientación de la elipse de error a partir de los autovalores y autovectores.

% 4. Elipse de Error Estándar
fprintf('--- ELIPSES DE ERROR ---\n');

% Autovalores y autovectores de la matriz de covarianza
[eig_vec, eig_val] = eig(Sx);
eigenvalores = diag(eig_val);

% Ordenar para identificar semieje mayor (s2) y menor (s1)
[eig_val_sorted, orden] = sort(eigenvalores, 'descend');
eig_vec_sorted = eig_vec(:, orden);

s2 = sqrt(eig_val_sorted(1)); % Semieje mayor
s1 = sqrt(eig_val_sorted(2)); % Semieje menor

% Orientación de la elipse (en grados centesimales o gons)
angulo_rad = atan2(eig_vec_sorted(1,1), eig_vec_sorted(2,1));
orientacion = mod(angulo_rad * 200 / pi, 200); 

fprintf('Semiejes estándar:\n');
fprintf('s1 = %.4f\n', s1);
fprintf('s2 = %.4f\n', s2);
fprintf('orientación = %.4f\n\n', orientacion);

5. Elipse de Error para Probabilidad del 91%

Cálculo de la elipse expandida utilizando un factor de escala c basado en una probabilidad específica.

% 5. Elipse de Error asociada a probabilidad 0.91
probabilidad = 0.91;

% P = 1 - exp(-c^2 / 2) -> despejando el factor de escala 'c'
c = sqrt(-2 * log(1 - probabilidad));

p1 = c * s1;
p2 = c * s2;

fprintf('Elipse (Prob = 0.91):\n');
fprintf('c = %.4f\n', c);
fprintf('p1 = %.4f\n', p1);
fprintf('p2 = %.4f\n', p2);