Mínimos cuadrados no lineales, basado en problemas
Este ejemplo muestra cómo realizar el ajuste de curvas de mínimos cuadrados no lineales con Flujo de trabajo de optimización basada en problemas.
Modelo
La ecuación modelo para este problema es
donde , , y son los parámetros desconocidos, es la respuesta y es el tiempo. Este problema requiere datos para tiempos tdata
y mediciones de respuesta (ruidosas) ydata
. El objetivo es encontrar el mejor valor de y , es decir, aquellos valores que minimizan
Muestrear datos
Habitualmente, tiene datos para un problema. En este caso, genere datos ruidosos artificiales para el problema. Utilice A = [1,2]
y r = [-1,-3]
como valores subyacentes, y utilice 200 valores aleatorios de 0
a 3 como datos de tiempo. Represente los puntos de datos resultantes.
rng default % For reproducibility A = [1,2]; r = [-1,-3]; tdata = 3*rand(200,1); tdata = sort(tdata); % Increasing times for easier plotting noisedata = 0.05*randn(size(tdata)); % Artificial noise ydata = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata) + noisedata; plot(tdata,ydata,'r*') xlabel 't' ylabel 'Response'
Los datos son ruidosos. Por ello, es probable que la solución no coincida demasiado bien con los parámetros originales A
y r
.
Enfoque basado en problemas
Para encontrar los parámetros A
y r
que mejor se ajusten, defina primero variables de optimización con esos nombres.
A = optimvar('A',2); r = optimvar('r',2);
Cree una expresión para la función objetivo, que es la suma de los cuadrados que se desean minimizar.
fun = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata); obj = sum((fun - ydata).^2);
Cree un problema de optimización con la función objetivo obj
.
lsqproblem = optimproblem("Objective",obj);
Para el enfoque basado en problemas, especifique el punto inicial como una estructura, con los nombres de variable como campos de la estructura. Especifique el valor A = [1/2,3/2]
inicial y el valor r = [-1/2,-3/2]
inicial.
x0.A = [1/2,3/2]; x0.r = [-1/2,-3/2];
Revise la formulación del problema.
show(lsqproblem)
OptimizationProblem : Solve for: A, r minimize : sum(arg6) where: arg5 = extraParams{3}; arg6 = (((A(1) .* exp((r(1) .* extraParams{1}))) + (A(2) .* exp((r(2) .* extraParams{2})))) - arg5).^2; extraParams{1}: 0.0139 0.0357 0.0462 0.0955 0.1033 0.1071 0.1291 0.1385 0.1490 0.1619 0.1793 0.2276 0.2279 0.2345 0.2434 0.2515 0.2533 0.2894 0.2914 0.2926 : : extraParams{2}: 0.0139 0.0357 0.0462 0.0955 0.1033 0.1071 0.1291 0.1385 0.1490 0.1619 0.1793 0.2276 0.2279 0.2345 0.2434 0.2515 0.2533 0.2894 0.2914 0.2926 : : extraParams{3}: 2.9278 2.7513 2.7272 2.4199 2.3172 2.3961 2.2522 2.1974 2.1666 2.0944 1.9566 1.7989 1.7984 1.7540 1.8318 1.6745 1.6874 1.5526 1.5229 1.5680 : :
Solución basada en problemas
Resuelva el problema.
[sol,fval] = solve(lsqproblem,x0)
Solving problem using lsqnonlin. Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
sol = struct with fields:
A: [2x1 double]
r: [2x1 double]
fval = 0.4724
Represente la solución resultante y los datos originales.
figure responsedata = evaluate(fun,sol); plot(tdata,ydata,'r*',tdata,responsedata,'b-') legend('Original Data','Fitted Curve') xlabel 't' ylabel 'Response' title("Fitted Response")
La gráfica muestra que los datos ajustados coinciden bastante bien con los datos ruidosos originales.
Compruebe lo bien que coinciden los parámetros ajustados con los parámetros originales A = [1,2]
y r = [-1,-3]
.
disp(sol.A)
1.1615 1.8629
disp(sol.r)
-1.0882 -3.2256
Los parámetros ajustados difieren en un 15% con respecto a A
y en un 8% con respecto a r
.
Las funciones no compatibles requieren fcn2optimexpr
Si la función objetivo no está compuesta por funciones elementales, deberá convertir la función en una expresión de optimización utilizando fcn2optimexpr
. Consulte Convert Nonlinear Function to Optimization Expression. Para el presente ejemplo:
fun = @(A,r) A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata); response = fcn2optimexpr(fun,A,r); obj = sum((response - ydata).^2);
El resto de los pasos para resolver el problema son idénticos. La única diferencia está en la rutina de representación, donde llama a response
en lugar de a fun
:
responsedata = evaluate(response,sol);
Para ver la lista de funciones compatibles, consulte Supported Operations for Optimization Variables and Expressions.