fsolve not varying input variable

I'm trying to use fsolve to solve for the best fit parameters to a model. The model is generated algorithmically, by compressing and rescaling the signal and time vectors, rather than using an analytic function. If the model is called directly, its output is indeed dependant upon all six input parameters. However, when fsolve is used to try to fit the model, it does not alter the last parameter.

Due to the operations in the model (which include piecewise polynomials and integration for normalizing), it is not compatible with the legal operations of an OptimizationExpression, which excludes the use of 'solve'.

I've changed the model inputs from a parameter structure to a vector of variables, and debugged by checking the values internally in the model function and residual function. There doesn't appear to be an overwriting issue. And yet, fsolve will only alter the first 5 parameters, and leaves the last one unchanged.

The model with parameter structure input (func_ADLv3), model with parameter vector input (func_ADLv31), and other functions used with the evaluation are attached.

The call script (functesting_RTDsuite) runs all of the functions in the suite, and attempts to fit the model parameters to a previously generated output.

Its a lot of code, but I think I'm doing something wrong with the fsolve code, or for some reason fsolve is determining that the last parameter does not affect the output, but direct calls of the model clearly show a dependance upon it.

Thanks for any help, or recommendations for other methods.

%% functesting_RTDsuite
% Test RTD function suite held in directory with example values to confirm usage
% Collin Schmidt 2024/01/12
clear all; close all; clc;

tpts = 500;
t = linspace(0,5,tpts); %(sec) Relevant time vector

%% func_L2
% Using values matching Logistic Inlet from RS1 simulations
P.A=30; %(1/sec)logistic func severity
P.C=0.3; %(sec)logistic time delay

[R1.E,R1.Epp,NormE1,R1.F,R1.Fpp,NormF1] = func_L2(t,P); %Generate RTD1 with L2 model

figure (1)
subplot(2,2,1)
title('func L2')
yyaxis left
plot(t,R1.E,'.',t,ppval(R1.Epp,t))
xlabel('Time(sec)')
ylabel('Instantaneous RTD (E)')
yyaxis right
plot(t,R1.F,'.',t,ppval(R1.Fpp,t))
ylabel('Cumulative RTD (F)')
axis([t(1),t(end),0,1.1*max(R1.F)])
legend('E','Epp','F','Fpp')

%% func_ADL
% Model parameters
P.A = 10;  %Severity, Base Logistic Function (<100)
P.B = 0.5;  %Compression extent, Left side (0,1]
P.B1 = 0.3; %Fraction of Left side compressed (0,1]
P.C = 1;  %Delay, Base Logistic Function (>0)
P.D = 10; %Stretch extent, Right side (0,?]
P.D1 = 0.8; %Fraction of Right side stretched (0,1]

% [R2.E,R2.Epp,NormE2,R2.F,R2.Fpp,NormF2] = func_ADLv2(P,t); %Generate R2 with ADLv1 model
tend = 5; %(sec) End time for fitting model. Must cover full convolution range
[R2.E,R2.Epp,NormE2,R2.F,R2.Fpp,NormF2,twp2,tend2,P2] = func_ADLv3(P,tend); %Generate RTD1 with L2 model
t2 = linspace(0,tend2,tpts); %(sec) Time vector for evaluating R2

subplot(2,2,2)
title('func ADLv1')
yyaxis left
plot(twp2,R2.E,'.',t2,ppval(R2.Epp,t2),'-')
xlabel('Time(sec)')
ylabel('Instantaneous RTD (E)')
axis([0,tend2,0,1.1*max(R2.E)])
yyaxis right
plot(twp2,R2.F,'.',t2,ppval(R2.Fpp,t2),'-')
ylabel('Cumulative RTD (F)')
axis([0,tend2,0,1.1*max(R2.F)])
legend('E','Epp','F','Fpp')

%% Convolution for RTD
[R3.E,R3.Epp] = func_RTDconvolution(R1.Epp,R2.Epp,t2); %Convolute RTDS R1,R2. Internally normalized

subplot(2,2,3)
title('func RTDConvolution')
yyaxis left
plot(t,ppval(R1.Epp,t),t2,ppval(R2.Epp,t2),t2,R3.E,'.',t2,ppval(R3.Epp,t2),'-')
xlabel('Time(sec)')
ylabel('Instantaneous RTD (E)')
axis([0,tend2,0,1.1*max(R1.E)])
legend('R1.Epp','R2.Epp','R3.E','R3.Epp')

%% Convolution for Arbitrary Input
I1.F = 3*ppval(R1.Fpp,t2); %Using scaled L2 cumulative distribution for input
I1.Fpp = pchip(t2,I1.F); %fit input to pchip piecewise polynomial
[O1,O1pp] = func_convolution(I1.Fpp,R3.Epp,t2); %Convolute arbitrary input with R3
[I1.E,I1.Epp,NormI1] = func_RTD_E_from_F(I1.Fpp,t2); %Convert input (arbitrarily) to instantaneous RTD

subplot(2,2,4)
title('func Convolution with Arbitrary Input')
yyaxis left
plot(t2,ppval(I1.Fpp,t2),t2,ppval(R3.Epp,t2),t2,O1,'.',t2,ppval(O1pp,t2),'-')
xlabel('Time(sec)')
ylabel('Instantaneous RTD (E)')
axis([0,tend2,0,1.1*max(I1.F)])
yyaxis right
plot(t2,I1.E,'.',t2,ppval(I1.Epp,t2),'-')
ylabel('E for Input Function')
legend('R1.Epp','R2.Epp','R3.E','R3.Epp')

%% Fiting with ADL
% Adjusting parameters off from target ADL
x0(1) = 11;  %Severity, Base Logistic Function (<100)
x0(2) = 0.5;  %Compression extent, Left side (0,1]
x0(3) = 0.4; %Fraction of Left side compressed (0,1]
x0(4) = 1.2;  %Delay, Base Logistic Function (>0)
x0(5) = 8; %Stretch extent, Right side (0,?]
x0(6) = 0.6; %Fraction of Right side stretched (0,1]
% x0 = [P.A, P.B, P.B1, P.C, P.D, P.D1]; %initial parameter guess
options = optimoptions('fsolve','InitDamping',1e-2);
fun = @(x)func_ResidSectionADL(x,R2.Fpp,twp2)

fun = function_handle with value:
    @(x)func_ResidSectionADL(x,R2.Fpp,twp2)

% %% Using Problem Structure Generation
% problem.options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt);
% problem.objective = @fun;
% problem.x0 = x0;
% problem.solver = 'fsolve';
% % Solve the problem.
% x = fsolve(problem)


norm((fun(x0+[0 0 0 0 0 sqrt(eps)])-fun(x0))/sqrt(eps))

ans = 0

%xout=fsolve(fun,x0,options);

%Evaluate ADL with optimized parameters
P4.A=xout(1); P4.B=xout(2); P4.B1=xout(3);

Unrecognized function or variable 'xout'.

P4.C=xout(4); P4.D=xout(5); P4.D1=xout(6);
[R4.E,R4.Epp,R4.NormE4,R4.F,R4.Fpp,R4.NormF4,twp4,tend4,P4] = func_ADLv3(P4,twp2(end));

figure(2)
plot(twp2,R2.E,'.',twp4,R4.E,'-')
xlabel('Time (sec)')
ylabel('Signal')
legend('Target','Fitted Function')

function [resid] = func_ResidSectionADL(x,ypp,t)
%Calculate Residual between section and ADL model
%Inputs:
%P - parameter structure describing ADL function
%y - signal from response data, normalized piecewise polynomial fit
%t - time vector for signal data

% %Pack up parameters
% P.A=abs(x(1)); P.B=abs(x(2)); P.B1=abs(x(3));
% P.C=abs(x(4)); P.D=abs(x(5)); P.D1=abs(x(6));

%Evaluate ADL model
[E,Epp,NormE,F,Fpp,NormF,twp,tend,POut] = func_ADLv31(x,t(end));

%Calculate Residual
resid = abs(ppval(Fpp,t)-ppval(ypp,t));
figure(10)
title ('Fsolve Progress')
yyaxis left
plot(t,ppval(ypp,t),'.',t,ppval(Fpp,t))
xlabel('Time')
ylabel('F')
yyaxis right
plot(t,resid)
ylabel('Residual')
legend('Signal','Fit Function','Residual')
drawnow
end

function [EOut,EppOut,NormEOut,FOut,FppOut,NormFOut,twp,tendOut,POut] = func_ADLv3(P,tend)
% Generate Instantaneous (E) and Cumulative (F) Residence time distribution
% curves from supplied parameters using Asymmetrically Dispersed Logicstic
% Model, which is developed algorithmically, rather than analytically.
% A base function (2parameter logistic function) is altered by compressing
% amplitude and stretching along time axis
% Collin Schmidt 2022/06/30

% Ouputs
% EOut,FOut do not have evenly spaced time points due to algorithmic
% approach. Evaluate piecewise polynomials (EppOut,FppOut) after model
% generation for specific timepoint spacing. Do not evaluate pp structures
% beyond model duration (tOut). Discrete output of EOut,FOut correspond to
% time vector (twp).

% Revision Control
% v1 - original attempt
% v2 - removed one time vector, consolidating
% v3 - changing time input to just end time, returning parameters and end
% time (which may be adjusted from input values internally to satisfy RTD
% requirements)

%% Inputs
%Parameter Structure
% P.A  %Severity, Base Logistic Function (<100)
% P.B  %Compression extent, Left side (0,1]
% P.B1 %Fraction of Left side compressed (0,1]
% P.C  %Delay, Base Logistic Function (>0)
% P.D  %Stretch extent, Right side (0,?]
% P.D1 %Fraction of Right side stretched (0,1]

%Time vector inputs
% tend - overall experiment duration, needs to cover full convolution range
% in later evaluation using this model

%% User Input: set resolution
tpts_base = 100; %Resolution for base 2 parameter logistic function

%% Enforce Parameter Constraints
%Severity, Base Logistic Function (<100)
P.A=abs(P.A); %force positive value
if P.A>25
    P.A=25;
end
%Compression extent, Left side (0,1]
P.B = abs(P.B); %force positive value
if P.B>0.99 %force maximum value
    P.B=0.98;
end
%Delay, Base Logistic Function (>0)
P.C = abs(P.C); %force positive value
%Stretch extent, Right side (0,?]
P.D = abs(P.D);
%Fraction of Right side stretched (0,1]
P.D1 = abs(P.D1);
if P.D1>0.99 %force maximum value
    P.D1=0.98;
end

%% Features of Base Logistic Function
% Left Edge of Base Function
delta = 0.9999; %(-)Wavelet signal threshold for edges
tle = P.C-log(delta/(1-delta))/P.A; %time of left edge of wavepacket
% Right Edge of Base Function
tre = P.C+log(delta/(1-delta))/P.A; %time of right edge of wavepacket
% Width of Base Function
w = abs(2/P.A*log((1-delta)/delta)); %(timeunits)Wavelet full width

%% **Automation of tw
% tw = linspace(tle,tre+w/2*P.D1*P.D,100); %enforce resolution of wavepacket
% %extends end of time accounting for stretch parameter

%% Features of Dispersal
tlc = P.C-w/2*(1-P.B1); %location where left side compression ends
trs = P.C-P.D1*(w/2)+w/2; %time location where right side stretch begins

%% Conditional Enforcement for left hand side
% wavelet must land completely in positive time domain after compression of left hand side
if P.C<w/2 %Check if base logistic function penetrates to negative times
    if P.B1<(w/2-P.C)/(w/2)
        P.B1 = (w/2-P.C)/(w/2); %enforce minimum fraction of left side compression
    end
    Bmin = tle/(tle-tlc); %min value of compression parameter for left side to terminate after t=0
    if P.B<Bmin
        P.B=Bmin %enforce minimum compression to prevent discontinuity at t=0
    end
end

%% Evaluate Base Function
tbase = linspace(tle,tre,tpts_base); %(timeunits) Time vector for base logistic function
tstep_base = tbase(2)-tbase(1); %(timeunits) Time vector step size for appending
[E_L2,Epp_L2,NormE_L2,F_L2,Fpp_L2,NormF_L2] = func_L2(tbase,P); %Parameter values A,C used for base logistic function

%% Left side partial compression
% Find max time location of left side compression
[M,Itlc] = min(abs(tbase-tlc)); %Find index in time vector for end of left side compression
% Left Side (Compression)
t_left  =tbase(1:Itlc); %Sample time vector, left side compression region
t_left = P.B*(tbase(Itlc)-t_left)+t_left; %shift timepoints left of center some fraction of their distance from the center
% Left Side Pad (constrain values to start of time vector)
padpts = 20;
t_left = [linspace(-w/2,0.9*t_left(1),padpts),t_left]; %append with points to meet zero
Ewp = [zeros(1,padpts),E_L2]; %append E vector with appropriate pad points on left hand side

%% Right side partial stretch
[M,Itrs] = min(abs(tbase-trs)); %Find index in time vector for start of right side stretch
% Right Side (Stretch)
t_right = tbase(Itrs:end); %Sample time vector, right side stretch region
a = 1.5; %1.5 gives smoother transition %stretch exponent on time displacement
t_right = P.D*(t_right-tbase(Itrs)).^a+t_right; %Stretched Time vector to skew function to the right
% Pad over complete time domain
% if t_right(end)<t(end) %determine if stretched time vector reaches full experiment domain
if t_right(end)<tend %if stretched wavepacket is within tend, pad with zeros to tend
    t_right = [t_right,linspace(t_right(end)+tstep_base,tend,padpts)]; %additional time points to reach desired end time
    Ewp = [Ewp,zeros(1,padpts)]; %append E vector with appropriate pad points on right hand side
    tendOut=t_right(end); %(timeunits) Save end time for output
else %stretched wave packet exceeding input 'tend', need to update output
    tendOut=t_right(end); %(timeunits) Save new end time for output
end %if

%% Central Region (between tlc,trs) unchanged
t_center = tbase(Itlc+1:Itrs-1); %Sample time vector, central unchanged region

%% Piecewise Polynomial Fit to Asymetrically Dispersed Logistic Function
twp = [t_left, t_center, t_right]; %dispersed and padded wavepacket time vector
% E_ADL=spline(twp,Ewp); %Spline for Instantaneous RTD (padded on left and right to constrain edge of spline)
E_ADL=pchip(twp,Ewp); %Spline for Instantaneous RTD (padded on left and right to constrain edge of spline)
[EOut,EppOut,NormEOut] = func_Normalize_pp_E(E_ADL,twp); %Normalize
[FOut,FppOut,NormFOut] = func_RTD_F_from_E(EppOut,twp); %Generate cumulative RTD (internally normalized)

% figure(12) %Comparison of Base Logistic Function and Asymmetrically Dispersed Logistic Function
% title('Base and Asymmetrically Dispersed Logistic Function')
% yyaxis left
% plot(tbase,E_L2,'.',twp,ppval(EppOut,twp))
% title('Wavepacket Domain')
% xlabel('Time')
% ylabel('E')
% yyaxis right
% plot(tbase,F_L2,'.',twp,ppval(FppOut,twp))
% axis([0 tOut 0 1.1])
% xlabel('Time')
% ylabel('F')
% legend({'Base E','ADL E','Base F','ADL F'},'Location','NorthEast')

POut=P; %save adjusted parameters for output
end %function
function [EOut,EppOut,NormEOut,FOut,FppOut,NormFOut,twp,tendOut,POut] = func_ADLv31(x,tend)
% Generate Instantaneous (E) and Cumulative (F) Residence time distribution
% curves from supplied parameters using Asymmetrically Dispersed Logicstic
% Model, which is developed algorithmically, rather than analytically.
% A base function (2parameter logistic function) is altered by compressing
% amplitude and stretching along time axis
% Collin Schmidt 2022/06/30

% Ouputs
% EOut,FOut do not have evenly spaced time points due to algorithmic
% approach. Evaluate piecewise polynomials (EppOut,FppOut) after model
% generation for specific timepoint spacing. Do not evaluate pp structures
% beyond model duration (tOut). Discrete output of EOut,FOut correspond to
% time vector (twp).

% Revision Control
% v1 - original attempt
% v2 - removed one time vector, consolidating
% v3 - changing time input to just end time, returning parameters and end
% time (which may be adjusted from input values internally to satisfy RTD
% requirements)
% v31 - changing input parameter P to vector x in an attempt to debug
% fsolve not iterating P.D1

%% Inputs
%Parameter Structure
% x(1)  %Severity, Base Logistic Function (<100)
% x(2)  %Compression extent, Left side (0,1]
% x(3) %Fraction of Left side compressed (0,1]
% x(4)  %Delay, Base Logistic Function (>0)
% x(5)  %Stretch extent, Right side (0,?]
% x(6) %Fraction of Right side stretched (0,1]

%Time vector inputs
% tend - overall experiment duration, needs to cover full convolution range
% in later evaluation using this model

%% User Input: set resolution
tpts_base = 100; %Resolution for base 2 parameter logistic function

%% Enforce Parameter Constraints
%Severity, Base Logistic Function (<100)
x(1)=abs(x(1)); %force positive value
if x(1)>25
    x(1)=25;
end
%Compression extent, Left side (0,1]
x(2) = abs(x(2)); %force positive value
if x(2)>0.99 %force maximum value
    x(2)=0.98;
end
%Delay, Base Logistic Function (>0)
x(4) = abs(x(4)); %force positive value
%Stretch extent, Right side (0,?]
x(5) = abs(x(5));
%Fraction of Right side stretched (0,1]
x(6) = abs(x(6));
if x(6)>0.99 %force maximum value
    x(6)=0.98;
end

%% Features of Base Logistic Function
% Left Edge of Base Function
delta = 0.9999; %(-)Wavelet signal threshold for edges
tle = x(4)-log(delta/(1-delta))/x(1); %time of left edge of wavepacket
% Right Edge of Base Function
tre = x(4)+log(delta/(1-delta))/x(1); %time of right edge of wavepacket
% Width of Base Function
w = abs(2/x(1)*log((1-delta)/delta)); %(timeunits)Wavelet full width

%% **Automation of tw
% tw = linspace(tle,tre+w/2*x(6)*x(5),100); %enforce resolution of wavepacket
% %extends end of time accounting for stretch parameter

%% Features of Dispersal
tlc = x(4)-w/2*(1-x(3)); %location where left side compression ends
trs = x(4)-x(6)*(w/2)+w/2; %time location where right side stretch begins

%% Conditional Enforcement for left hand side
% wavelet must land completely in positive time domain after compression of left hand side
if x(4)<w/2 %Check if base logistic function penetrates to negative times
    if x(3)<(w/2-x(4))/(w/2)
        x(3) = (w/2-x(4))/(w/2); %enforce minimum fraction of left side compression
    end
    Bmin = tle/(tle-tlc); %min value of compression parameter for left side to terminate after t=0
    if x(2)<Bmin
        x(2)=Bmin %enforce minimum compression to prevent discontinuity at t=0
    end
end

%% Evaluate Base Function
tbase = linspace(tle,tre,tpts_base); %(timeunits) Time vector for base logistic function
tstep_base = tbase(2)-tbase(1); %(timeunits) Time vector step size for appending
P.A=x(1); P.C=x(4); %L2 parameters from input vector
[E_L2,Epp_L2,NormE_L2,F_L2,Fpp_L2,NormF_L2] = func_L2(tbase,P); %Parameter values A,C used for base logistic function

%% Left side partial compression
% Find max time location of left side compression
[M,Itlc] = min(abs(tbase-tlc)); %Find index in time vector for end of left side compression
% Left Side (Compression)
t_left  =tbase(1:Itlc); %Sample time vector, left side compression region
t_left = x(2)*(tbase(Itlc)-t_left)+t_left; %shift timepoints left of center some fraction of their distance from the center
% Left Side Pad (constrain values to start of time vector)
padpts = 20;
t_left = [linspace(-w/2,0.9*t_left(1),padpts),t_left]; %append with points to meet zero
Ewp = [zeros(1,padpts),E_L2]; %append E vector with appropriate pad points on left hand side

%% Right side partial stretch
[M,Itrs] = min(abs(tbase-trs)); %Find index in time vector for start of right side stretch
% Right Side (Stretch)
t_right = tbase(Itrs:end); %Sample time vector, right side stretch region
a = 1.5; %1.5 gives smoother transition %stretch exponent on time displacement
t_right = x(5)*(t_right-tbase(Itrs)).^a+t_right; %Stretched Time vector to skew function to the right
% Pad over complete time domain
% if t_right(end)<t(end) %determine if stretched time vector reaches full experiment domain
if t_right(end)<tend %if stretched wavepacket is within tend, pad with zeros to tend
    t_right = [t_right,linspace(t_right(end)+tstep_base,tend,padpts)]; %additional time points to reach desired end time
    Ewp = [Ewp,zeros(1,padpts)]; %append E vector with appropriate pad points on right hand side
    tendOut=t_right(end); %(timeunits) Save end time for output
else %stretched wave packet exceeding input 'tend', need to update output
    tendOut=t_right(end); %(timeunits) Save new end time for output
end %if

%% Central Region (between tlc,trs) unchanged
t_center = tbase(Itlc+1:Itrs-1); %Sample time vector, central unchanged region

%% Piecewise Polynomial Fit to Asymetrically Dispersed Logistic Function
twp = [t_left, t_center, t_right]; %dispersed and padded wavepacket time vector
% E_ADL=spline(twp,Ewp); %Spline for Instantaneous RTD (padded on left and right to constrain edge of spline)
E_ADL=pchip(twp,Ewp); %Spline for Instantaneous RTD (padded on left and right to constrain edge of spline)
[EOut,EppOut,NormEOut] = func_Normalize_pp_E(E_ADL,twp); %Normalize
[FOut,FppOut,NormFOut] = func_RTD_F_from_E(EppOut,twp); %Generate cumulative RTD (internally normalized)

% figure(12) %Comparison of Base Logistic Function and Asymmetrically Dispersed Logistic Function
% title('Base and Asymmetrically Dispersed Logistic Function')
% yyaxis left
% plot(tbase,E_L2,'.',twp,ppval(EppOut,twp))
% title('Wavepacket Domain')
% xlabel('Time')
% ylabel('E')
% yyaxis right
% plot(tbase,F_L2,'.',twp,ppval(FppOut,twp))
% axis([0 tOut 0 1.1])
% xlabel('Time')
% ylabel('F')
% legend({'Base E','ADL E','Base F','ADL F'},'Location','NorthEast')

% POut=P; %save adjusted parameters for output
POut.A=x(1);
POut.B=x(2);
POut.B1=x(3);
POut.C=x(4);
POut.D=x(5);
POut.D1=x(6);

end %function
function [O,Opp] = func_convolution(Ipp,RTDpp,t)
% Return convolution of arbitrary Input (Ipp, piecewise polynomial) with
% normalized instantaneous RTD (RTDpp, piecewise polynomial). **Output is
% in units of input function, and is NOT normalized.
% Discreet vector of time (t) must span full convolution.

O = zeros(1,length(t)); %Instantiate
for it=1:length(t)
    T=t(it); %(sec)Time to evaluate outlet flow
    %Convolution integral to calculate output
    %Options: Ignore relative tolerance, use absolute tolerance value
    O(it) = integral(@(tau) (ppval(Ipp,tau).*ppval(RTDpp,T-tau)),0,T,'RelTol',0,'AbsTol',1e-6);
end
% O(isnan(O))=0; %Set NaN values in 'O' to zero
Opp = spline(t,O); %Piecewise polynomial fit to data
end %function
function [EOut,EppOut,NormE,FOut,FppOut,NormF] = func_L2(t,P)
%Returns instantaneous and cumulative RTDs for 2 parameter logistic model.

%Inputs
% P.A - (1/timeunits) logistic function severity
% P.C - (timeunits) delay
% t - (timeunits) must be linearly spaced for gradient calculation of
%instantaneous RTD.

F = 1./(1+exp(-P.A*(t-P.C))); %analytic evaluation
Fpp = pchip(t,F); %piecewise polynomial fit for cumulative RTD (pchip)
[FOut,FppOut,NormF] = func_Normalize_pp_F(Fpp,t); %normalize discrete and pp of cumulative RTD
[EOut,EppOut,NormE] = func_RTD_E_from_F(FppOut,t); %convert to instantaneous RTD (normalizes internally)
end %function
function [EOut,EppOut,NormE] = func_Normalize_pp_E(Epp,t)
%Normalize discreet and piecewise polynomial structure for instantaneous RTD
NormE = integral(@(tau) ppval(Epp,tau),t(1),t(end),'RelTol',0,'AbsTol',1e-6); %current norm integral
EOut=ppval(Epp,t)/NormE; %Normalize
EppOut = spline(t,EOut); %fit to normalized instantaneous RTD
NormE = integral(@(tau) ppval(EppOut,tau),t(1),t(end),'RelTol',0,'AbsTol',1e-6); %normalization factor for output
end %function
function [FOut,FppOut,NormF] = func_Normalize_pp_F(Fpp,t)
%Normalize discreet and piecewise polynomial strucuture for cumulative RTD or Step Response
NormF = ppval(Fpp,t(end)); %take steady state value as normalization
FOut=ppval(Fpp,t)/NormF; %Normalize
FppOut = pchip(t,FOut); %fit to normalized instantaneous RTD
NormF = ppval(FppOut,t(end)); %check normalization
end %function
function [EOut,EppOut,NormE] = func_RTD_E_from_F(Fpp,t)
%Instantaneous RTD from Cumulative RTD by taking gradient and normalizing
%Time vector requires regular linear spacing
E = gradient(ppval(Fpp,t));
Epp = spline(t,E); %Piecewise polynomial fit to raw instantaneous RTD
[EOut,EppOut,NormE] = func_Normalize_pp_E(Epp,t); %Normalization on Instantaneous RTD
NormE; %normalization check. Must reach a max of 1.
end %function
function [FOut,FppOut,NormF] = func_RTD_F_from_E(Epp,t)
%Cumulative RTD from Instantaneous RTD by integration and normalizing
F = zeros(1,length(t)); %Instantiate
for it=1:length(t)
    T=t(it); %(sec)Time to evaluate outlet flow
    F(it) = integral(@(tau) ppval(Epp,tau),t(1),T,'RelTol',0,'AbsTol',1e-6);
end
Fpp = pchip(t,F); %Piecewise polynomial fit to cumulative RTD
[FOut,FppOut,NormF] = func_Normalize_pp_F(Fpp,t);
NormF; %normalization check. Must reach a max of 1.
end %function
function [EOut,EppOut] = func_RTDconvolution(RTDpp_1,RTDpp_2,t)
% Return normalized Instantaneous RTD for convolution of two instantaneous
% RTDs. **Both input RTDs must be normalized, and fitted over full time
% range.

O = zeros(1,length(t)); %Instantiate 
for it=1:length(t)
    T=t(it); %(sec)Time to evaluate outlet flow
    %Convolution integral to calculate output RTD
    %Options: Ignore relative tolerance, use absolute tolerance value
    O(it) = integral(@(tau) (ppval(RTDpp_1,tau).*ppval(RTDpp_2,T-tau)),0,T,'RelTol',0,'AbsTol',1e-6);
end
% O(isnan(O))=0; %Set NaN values in 'O' to zero
Opp = spline(t,O); %Piecewise polynomial fit to resulting instantaneous RTD
[EOut,EppOut,NormE] = func_Normalize_pp_E(Opp,t); %Normalize fitted instantaneous RTD
NormE; %Check normalization
end %function

Answer 1

[Torsten]

I included your code and inserted the line

norm((fun(x0+[0 0 0 0 0 sqrt(eps)])-fun(x0)))

before your call to "fsolve".

sqrt(eps) is the default FiniteDifferenceStepSize to compute the sensitivity of your system of equations with respect to the parameters. I get 0 as answer. This is bad because if the sensitivity is 0, "fsolve" won't vary the parameter because it gets signaled that your model does not depend on it.

So give ascending values to sqrt(eps) in the expression, see what MATLAB returns and set

FiniteDifferenceStepSize

in the optimoptions structure to the value where you first get a sensitivity different from 0.

As @Matt J describes below, such lack of sensitivity often is an indication that your functions do not depend differentiably on the parameters. This is often caused by if-statements and other artificial interventions in the computation of the function values.

[Collin Schmidt]

This worked. Also used suggestions from other commentors. In summary:

1) Switched to lsqnonlin, better match to programming problem than fsolve

2) confirmed from lsqnonlin output 'jacobian' that step in last parameter had zero impact on residual

3) increased FiniteDifferenceStepSize for last parameter manually until Jacobian was non-zero, and parameter was iterated. Other terms in step size vector were set as "sqrt(eps)", was unsure how else to keep them at default values.

Final code attached.

[Matt J]

increased FiniteDifferenceStepSize for last parameter manually until Jacobian was non-zero

Or maybe the final parameter just needs to be scaled, i.e., expressed in different units, so that the residuals are more sensitive to it.

Answer 2

[Matt J]

I can't say definitively whether this is responsible for the specific problem you cite, but you are not being at all careful in your implementation to avoid discontinuities and non-differentiabilities, thus violating theoretical assumptions of fsolve.

This for example,

if x(2)>0.99 %force maximum value
    x(2)=0.98;
end

will introduce a jump discontinuity.

This,

x(5) = abs(x(5));

will not introduce a discontinuity, and so is not as impactful, but it does introduce non-differentiabilityat x(5)=0.

You are also creating potential hazards with things like this,

if x(1)>25
    x(1)=25;
end

which creates a very large plateau region where the derivatives with respect to x(1) are 0. If the solver lands there, the iterations could get stuck (though it depends on how strongly x(1) is coupled with the other parameters).

What you really should be doing is using a different solver which allows you to impose bounds, e.g., lsqnonlin, and maybe more complicated constraints.

[Collin Schmidt]

That's very helpful, and I think it is pointing towards the problem.

An analytic function is used to generate a logistic function for supplied time values. The model (ADLv31) uses the last parameter to rescale a time vector, effectively stretching the function over a given range.

The reshaped discrete function is used in the evaluation of the residual, so I assumed that fsolve would be able to determine the derivative with respect to each parameter by tweaking the values and looking at the residual. From what you are saying, I am incorrect in that, and I'm assuming that the linear operations on the time vector result in the derivative being uncalculatable. Do I have that right?

[Collin Schmidt]

I implemented lsqnonlin with the parameter constraints as you suggested (fit_lsqnonlin_ADLv33). The abs value actions and parameter forcing has been removed from the model (func_ADLv33).

However, the last parameter is still not being altered by the solver.

If I feed the correct value for the last parameter (data to fit is from a call of func_ADLv33, so the correct parameter values are known), the optimization proceeds correctly, and accurately fits the rest of the parameters.

If I do not feed the correct value for the last parameter, lsqnonlin is not able to fit the data well, and monitoring the progress shows all the parameters but the last one being altered.

Does that help to narrow it down at all?

[Torsten]

One way is to choose a higher value for the "FiniteDifferenceStepSize" in the options setting for the 6th parameter.

A better way is to find the reason in your code why a too small change in the 6th fit parameter does not lead to a change in the values returned by "func_ResidSectionADL". If you find and you can remedy this reason, the parameter will change in the course of the solution process.

I think you know best where decimal places are wiped out in your code.

[Collin Schmidt]

Other than the parameter forcing (which was removed when using lsqnonlin), I do not believe there are places where decimals are wiped out. There are some numeric integration routines, but they are executed at higher resolution. Also, (qualitatively) changing the last parameter by even 10% results in a major change to the function, so its odd to me what the sensitivity is coming back as zero.

However, when I changed the solver to 'fminsearch', all parameters are altered and the solver is able to converge. I tried fminsearch because it is a derivative free method, following the comments above that fsolve and lsqnonlin would not be able to interpret derivatives from the model as written.

fminsearch is an inferior method, and the fitting is not excellent, so I'm going to try to supply the Jacobian with the "SpecifyObjectiveGradient" option on lsqnonlin.

[Torsten]

Also, (qualitatively) changing the last parameter by even 10% results in a major change to the function, so its odd to me what the sensitivity is coming back as zero.

The sensitivity is coming back as zero because the change of the parameter by sqrt(eps) does not cause your function to return something different compared to the call with the parameter unchanged.

Nonetheless, you have to find the reason for this "oddity" because it is the reason why your parameter is not changed. Or you have to adapt the sensitivity of this parameter by changing the "FiniteDifferenceStepSize" for this parameter in the options setting. All methods except for "fminsearch" rely on the Jacobian matrix, and the respective solvers get signalled that the entry in the Jacobian is zero if the values returned by "func_ResidSectionADL" do not change. By the way: the parameter in "fsolve" is also altered, namely by sqrt(eps), but since the feedback remains unchanged, also the parameter is not changed in the following Newton-step.

fminsearch is an inferior method, and the fitting is not excellent, so I'm going to try to supply the Jacobian with the "SpecifyObjectiveGradient" option on lsqnonlin.

Why do you want to specify the Jacobian ? You will have to use a finite-difference approximation as lsqnonlin does itself. Do you think you can do better ?

[Matt J]

Other than the parameter forcing (which was removed when using lsqnonlin), I do not believe there are places where decimals are wiped out.

What is the Jacobian and the final residual vector? I recommend you run lsqnonlin with all of its output arguments

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(___)

and post htem in a .mat file.