Wielu badaczy, jeszcze w czasach gdy obliczenia były zaporowo drogie, zaproponowało rozmaite transformacje powyższego problemu, pozwalające w ostatecznym rozrachunku sformułować zadanie jako liniowe zadanie najmniejszych kwadratów. Mając w pamięci numeryczną maksymę, że

Jeśli zadanie jest za trudne, należy… zmienić zadanie!

możemy ulec wrażeniu, że będzie to właściwa droga ku efektywnemu znalezieniu ,,optymalnych” wartości szukanych parametrów.

W przypadku zadania dopasowania Km i V, używano (por. [10]) między innymi takich (matematycznie równoważnych) transformacji (4.1):

Solversqp jest narzędziem służącym do rozwiązywania znacznie bardziej złożonych zadań optymalizacji, dlatego jego użycie bywa skomplikowane. Ale w naszym prostym przypadku — minimalizacji bez ograniczneń, dla funkcjonału kwadratowego — nie będziemy musieli podawać mu zbyt wielu parametrów. Najpierw jednak zdefiniujemy funkcjonał, który będziemy minimalizować:

% (X(1), X(2)) <---> (V, K)
F = @(X,S) (X(1)*S)./(X(2)+S);
phi = @(X) sumsq(v0 - F(X,S));

Jak widać, zaczęliśmy od określenia funkcji

która pojawia się w definicji (4.1) i w (4.3), i która z pewnością nam się przyda: na przykład do narysowania wykresu F dla wyznaczonych V i Km (dlatego definiujemy ją od razu wektorowo ze względu na S). Potem, przy użyciu F, zdefiniowaliśmy ϕ, przy czym — ponieważ phi jest anonimową funkcją tylko jednego argumentu, X, to w jej definicji zostaną użyte wcześniej przez nas zdefiniowane wektory S i v0, zawierające dane zadania. Dzięki temu, łatwo nam będzie wyznaczyć sumę kwadratów współrzędnych wektora v0 - F(X,S) — do tego właśnie służy funkcja Octave sumsq.

Teraz wystarczy wywołać solver sqp na naszym funkcjonale ϕ, z początkowym przybliżeniem V,KmT wyznaczonym np. metodą Dixona:

% podane na wejściu sqp parametry (V, Km) zostały wyznaczone metodą Dixona
[x, phix, info, iter, nphi] = sqp([V;Km], phi)
V = x(1); Km = x(2);

W wyniku dostajemy V=25.40 oraz Km=2.33, dla których ϕ⁢V,Km=1.51, a więc (nieco) lepiej niż dotychczas! Jak zwykle w przypadku metod iteracyjnych, zażądamy możliwie wielu informacji diagnostycznych, takich jak: kod zakończenia info, liczba wykonanych iteracji iter, oraz wywołań funkcjonału nphi. Pożądana wartość info to, zgodnie z manualem, 101: zakończenie z powodu nikłego postępu iteracji. Musimy jednak pamiętać, że znalezione minimum może być jedynie minimum lokalnym i wybór innego punktu startowego może dać w rezultacie znacznie lepszy (lub gorszy) wynik.

Podsumowując, w dzisiejszych czasach, gdy nawet dość zaawansowae obliczenia są (w miarę) łatwe i tanie, nie powinniśmy bać się nieliniowości. Pamiętajmy także, że jeśli tylko możemy sensownie wspomóc się przybliżeniem uzyskanym na drodze rozsądnej linearyzacji — warto z tej opcji skorzystać.

Na marginesie dodajmy, że opisane przez nas wcześniej metody sprowadzenia zadania do zadania liniowego, są wciąż rutynowo stosowane w innych zadaniach dopasowywania do punktów pomiarowych wykresu, na przykład postaci y⁢x=a⁢eb⁢x: biorąc logarytm od obu stron, ponownie dostajemy zadanie liniowe na a i b: log⁡y=log⁡a+b⁢x. Jednak, jak już doświadczyliśmy, transformacja zadania najczęściej powoduje też zmianę sposobu, w jaki mierzymy błąd dopasowania. Niektóre transformacje dodatkowo mogą wyolbrzymiać lub redukować różnice pomiędzy danymi punktami pomiarowymi.

S = [0.25 ; 0.30; 0.40; 0.50; 0.70; 1.00; 1.40; 2.00];
v0 = [2.4; 2.6; 4.2; 3.8; 6.2; 7.4; 10.2; 11.4];

F = @(X,S) (X(1)*S)./(X(2)+S);
fitfun = @(X) sumsq(v0 - F(X,S));

function cbA = lsqcond(A, b, x, tol)
if nargin < 4
	tol = 1e2;
end
c = norm(A*x)/norm(b);
t = norm(A*x-b)/norm(A*x);
condA = cond(A);
cb = condA/c;
cA = (condA^2)*t + condA;
cbA = cb+cA;
if (cbA) > tol
	warning(['Uwarunkowanie zadania: ', num2str([cb, cA])]);
end
end

p = 0; n = 100;
s = linspace(0.125,4,n);

e = ones(size(S));

model = 5;
A = [1./S, -1./v0]; b = -e;
x = A \ b;
V(model) = x(2); Km(model) = x(1);
resid(model) = fitfun([V(model),Km(model)]);
cbA(model) = lsqcond(A,b,x);
plot(A(:,1), A(:,2), 'ro'); hold on;
plot(1./s, -(1+x(1)./s)/x(2),'-b'); hold off;
title('Fit: liniowy'); pause(p);
print('-depsc','linear.eps');

% Lineweaver i Burk
model = 1;
A = [1./S, e]; b = 1./v0;
x = A \ b;
V(model) = 1/x(2); Km(model) = x(1)*V(model);
resid(model) = fitfun([V(model),Km(model)]);
cbA(model) = lsqcond(A,b,x);
plot(A(:,1), b, 'ro'); hold on;
plot(1./s, x(1)./s+x(2),'-b'); hold off;
title('Fit: Lineweaver i Burk'); pause(p);
print('-depsc','lib.eps');

% Dixon
model = 2;
A = [e, S]; b = S./v0;
x = A \ b;
V(model) = 1/x(2); Km(model) = x(1)*V(model); 
resid(model) = fitfun([V(model),Km(model)]);
cbA(model) = lsqcond(A,b,x);
plot(A(:,2), b, 'ro'); hold on;
plot(s, x(2)*s+x(1),'-b'); hold off;
title('Fit: Dixon'); pause(p);
print('-depsc','dixon.eps');


% Eadie i Hofstee
model = 3;
A = [-v0./S, e]; b = v0;
x = A \ b;
V(model) = x(2); Km(model) = x(1);
resid(model) = fitfun([V(model),Km(model)]);
cbA(model) = lsqcond(A,b,x);
plot(A(:,1), b, 'ro'); hold on;
s3 = linspace(min(A(:,1))-2, max(A(:,1))+2, 100);
plot(s3, x(1)*s3+x(2),'-b'); hold off;
title('Fit: Eadie i Hofstee'); pause(p);
print('-depsc','eih.eps');

% nonlinear fit; X = (V,Km)
model = 4;
[x, fitfunx, info, iter, nF, lambda] = sqp([V(2);Km(2)], fitfun);
V(model) = x(1); Km(model) = x(2); 
resid(model) = fitfun([V(model),Km(model)]);
cbA(model) = NA;
if info ~= 101
	warning('SQP: problemy z rozwiazaniem zadania');
end

printf('+------+---------+--------+----------+----------+\n');
printf(' Model | V       | K      | resid    | cond\n');
printf('+------+---------+--------+----------+----------+\n');
printf(' %2d    | %4.2f   | %4.2f   | %6.2f   | %4.2e\n', [1:length(V);V;Km;resid;cbA]);
printf('+------+---------+--------+----------+----------+\n');

plotname = {'dane','Lineweaver i Burk','Dixon','Eadie i Hofstee','NLZNK','liniowy'};
plot(S, v0, 'ro');
hold on;
for i = 1:length(V)
	plot(s, F([V(i),Km(i)],s) , ['-', num2str(i-1)]);
	pause(p);
end
hold off;
xlabel('S'); ylabel('v_0');
legend(plotname, 'location', 'southeast');
print('-depsc','michaelis-fit.eps')