enforce smooth parameter variation in model fit across space

Genetic algorithm ga() was used to fit a model to data taken across 3D space.

In this problem, 6 parameters are fit for each spatial point, which takes up to 2 minutes on a laptop; there are 2000 points. Fitting the model to each and every point sequentially takes 1-2 days.

As a first-pass, each point was considered independently, even though the physics-based model parameters must vary smoothly over space. What is a computationally-efficient way to force this smooth spatial variation on all 6 parameters (ideally without requiring a supercomputer or several days of processing)?

The 6 parameters only vary by about one order of magnitude, but I'm not able to say that their variation across XYZ space can be described by some simple function. At present, an optimal fit of the 6 parameters is obtained using ga() to minimize an objective function, as below:

function optimalParamsfor1Point = DoGeneticAlgFit(RawDatafor1Point,BackgroundData)
% setup steps ...
% ----------- DO THE OPTIMIZATION -----------------------------------------
options = optimoptions('...');
optimalParamsfor1Point = ga(@objfun,6,A,b,Aeq,beq,lb,ub,nlcon,options);
% -------------------------------------------------------------------------
    function cost = objfun(x)
        % given 6 element vector x, RawDatafor1Point, and BackgroundData, return a cost.
        % Low cost means parameters in vector x accurately describe raw data.
        ...
            cost = RMS(rawdatvstime - modelestimatevstime)
    end
end

So existing code has this structure (N=2000):

% inputs:
% RawXYZData is a Nx1 cell of recorded data vs time for each point
% XYZ is a Nx3 matrix of known XYZ coordinates
% BackgroundData is a struct containing fixed constants
% output:
optimalParams = nan(N,6);   % matrix will hold best fit parameters for each point
for i = 1:N
    optimalParams(XYZpoint,:) = DoGeneticAlgFit(RawXYZData{i},BackgroundData)
end

A solution to the whole problem has Nx6 (12,000) parameters. I have a solution now composed of N individual solutions, but I would like to impose another constraint: smooth variation across XYZ space.

More description on what I am doing:

I have used a 3D mapper system to measure interactions between a set of antennas. The mapper moves to a point and plays an excitation signal for some time while responses are measured across the system. So, at each point I have a time series of excitation and responses.

An electrical circuit model, with some unknown parameters that depend on the position, describes the interactions. Knowing the model, the excitation, and the responses I can identify best fit parameters.

Seems that I should just look at the shape of parameter solutions for independently-solved points and then pick a function (example, polynomial in X, Y, & Z) and solve globally for the terms of this function. To me this appears to be a messy undertaking: Suppose I assume a quartic equation in 3 space variables (35 terms) and then solve for those terms across my 6 parameters. This would be optimizing 210 unknowns simultaneously.

Is there a more straightforward approach?

Answer 1

[Star Strider]

It likely depends on your problem. If it involves integrating a set of nonlinear ordinary differential equations, and if the parameters span a large range of orders-of-magnitude, use a stiff solver such as ode15s or ode23s instead of ode45.

Does your problem optimise the parameters, or fit them to data?

I would consider the entire model together, and estimate all the parameters simultaneously.

[Brett]

Thanks Star Strider!

The 6 parameters only vary by about one order of magnitude, but I'm not able to say that their variation across XYZ space can be described by some simple function. At present, an optimal fit of the 6 parameters is obtained using ga() to minimize an objective function, as below:

function optimalParamsfor1Point = DoGeneticAlgFit(RawDatafor1Point,BackgroundData)
    % setup steps ...
    % ----------- DO THE OPTIMIZATION -----------------------------------------
    options = optimoptions('...');
    optimalParamsfor1Point = ga(@objfun,6,A,b,Aeq,beq,lb,ub,nlcon,options);
    % -------------------------------------------------------------------------
    function cost = objfun(x)
        % given 6 element vector x, RawDatafor1Point, and BackgroundData, return a cost.
        % Low cost means parameters in vector x accurately describe raw data.
        ...
        cost = RMS(rawdatvstime - modelestimatevstime)
    end
end

So existing code has this structure (N=2000):

% inputs:
% RawXYZData is a Nx1 cell of recorded data vs time for each point
% XYZ is a Nx3 matrix of known XYZ coordinates
% BackgroundData is a struct containing fixed constants

% output:
optimalParams = nan(N,6);   % matrix will hold best fit parameters for each point
for i = 1:N
    optimalParams(XYZpoint,:) = DoGeneticAlgFit(RawXYZData{i},BackgroundData)
end

A solution to the whole problem has Nx6 (12,000) parameters. I have a solution now composed of N individual solutions, but I would like to impose another constraint: smooth variation across XYZ space.

Its not clear to me how this would be accomplished using a solver like ode45. Could you elaborate?

[Torsten]

What kind of model do you use to compute "modelestimatevstime" ? An algebraic relation ? A differential equation ?

[Star Strider]

I'm not certain that I understand what you're doing.

I've done similar optimisations, one being Parameter Estimation for a System of Differential Equations

The essence of that is this code --

t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

c=[0.902	0.06997	0.02463	0.00218
0.8072	0.1353	0.0482	0.008192
0.6757	0.2123	0.0864	0.0289
0.5569	0.2789	0.1063	0.06233
0.4297	0.3292	0.1476	0.09756
0.3774	0.3457	0.1485	0.1255
0.2149	0.3486	0.1821	0.2526
0.141	0.3254	0.194	0.3401
0.04921	0.2445	0.1742	0.5277
0.0178	0.1728	0.1732	0.6323
0.006431	0.1091	0.1137	0.7702
0.002595	0.08301	0.08224	0.835]; 
% theta0=[1;1;1;1;1;1];

% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 200;
Parms = 10;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);   % With 'HybridFcn'

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2025-12-18 23:23:14.2755

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2025-12-18 23:25:40.2996

GA_Time = etime(t1,t0)

GA_Time = 146.0241

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

Elapsed Time:   1.460241240000000E+02		00:02:26.024

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 0.1462
Generations 				 = 4257

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.32986
		Theta( 2) =  0.78067
		Theta( 3) =  4.74999
		Theta( 4) =  9.92945
		Theta( 5) =  0.83910
		Theta( 6) =  0.14558
		Theta( 7) =  0.94421
		Theta( 8) =  0.07683
		Theta( 9) =  0.00700
		Theta(10) =  0.00001

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dC=dcdt;
    end
C=Cv;
end

This will produce the best results if 'PopSz' is larger, however it will take longer to converge.

.

[Brett]

To Torsten's question:

"What kind of model do you use to compute "modelestimatevstime" ?"

It is an algebraic relation, but not a simple closed form equation.

Thank you Star Strider for the example code:

Your solution is very close to what I already have at each point in space. Now I want an efficient way to merge all points so that my parameters (your "Theta") vary smoothly between nearby points.

More description on what I am doing:

I have used a 3D mapper system to measure interactions between a set of antennas. The mapper moves to a point and plays an excitation signal for some time while responses are measured across the system. So, at each point I have a time series of excitation and responses.

An electrical circuit model, with some unknown parameters that depend on the position, describes the interactions. Knowing the model, the excitation, and the responses I can identify best fit parameters.

Seems that I should just look at the shape of parameter solutions for independently-solved points and then pick a function (example, polynomial in X, Y, & Z) and solve globally for the terms of this function. To me this appears to be a messy undertaking: Suppose I assume a quartic equation in 3 space variables (35 terms) and then solve for those terms across my 6 parameters. This would be optimizing 210 unknowns simultaneously.

Is there a more straightforward approach?

[Star Strider]

My pleasure!

I 'sort of' understand what you're doing (I'm an amateur radio operator). I don't entirely because I don't understand how you're defining the points for the exciter or antenna array. (In my chemical kinetics model, the data are functions only of time, and the model fits the equations by varying the parameter values, integrating the equations, and comparing the integrated equation results to the data.) I'm not certain that I understand how you're defining the parameters in your spatial model, although I suspect some sort of inverse-square approach might be necesary. Are you determining all the parameters for the entire antenna array for each exciter position? If so, that may be the best you can hope for. I'm not certain there would be a way to combine those estimates, or plot them with respect to the different exciter positions. If the antennas are omnidirectionsl, this task might be easier than if they were directional.

It would help to have some sort of description (in code as well as symbolic), and the appropriate data to fit it.