Developer Forum

ello,

I'm using simulated annealing (SA) to estimate parameters in a high-dimensional space (10--20 p). The code that I will show you (below) is a simplification of my programs but is useful for describing my problem.

SA requires that starting parameter estimates (seeds) be provided and bounds on these parameters are also useful to finding a solution. I have ten years of time series data for nonnative lake trout in Yellowstone Lake; in this case, the data are length-frequency histograms. The SA solver's objective is to minimize the least-squares fit of the observed (measured) length-frequency distribution to the distribution predicted by my model.

Because environmental conditions vary significantly among years, the seeds and their bounds differ among years. (Both types of data are preliminarily estimated using simple regression techniques.) Accordingly, I must send a year-specifc set of parameter seeds and their bounds to the solver with each time step (year). This is where my problem seems to come in.

Here's my MasterFunction code; I've put some of my addtional comments in quotes:

function [ParaEstsOut]=MasterFunc
randn('state',sum(100*clock));
rand('state',sum(100*clock));
% =========================================================================
fits = 3; "3 SA fits of the model to the data per year"
Paras = 6; "6 parameters to be estimated"
% =========================================================================
load ProjectYrs "load required data from current directory"
load LTGNMeshesMM
load LTBinTLIndex
load MethodsCount
load DataAnnLTLngFrqByMethodMesh
load DataAnnLTTLCountByMethodMesh
load DataQuadParameters "this is the file that contains data for parameter seeds and their bounds, indexed by year"
% =========================================================================
% Start with the default options
options = saoptimset;
% Modify options setting
options = saoptimset(options,'TemperatureFcn',@temperatureexp); % @temperatureexp @temperatureboltz @temperaturefast
options = saoptimset(options,'AnnealingFcn', @annealingfast); % @annealingboltz @annealingfast
options = saoptimset(options,'TolFun', 1e-8); % 1e-6
options = saoptimset(options,'StallIterLimit', 3000); % Default: 500 x number of parameters.
options = saoptimset(options,'MaxFunEvals', 180); % 3000*Paras
options = saoptimset(options,'Display', 'off');
options = saoptimset(options,'PlotFcns', { @saplottemperature @saplotbestx @saplotf @saplotstopping });
options = saoptimset(options,'OutputFcns', { [] });
options = saoptimset(options,'InitialTemperature', 250); % 100
options = saoptimset(options,'ReannealInterval', 50); % 100
% =========================================================================
ParaEstsOut=zeros((fits*(ProjectYrs-71)),Paras+5); % 1998 (72) is first year with data.
Counter = 1;
for yr = 72:ProjectYrs %:ProjectYrs % 1998 (72) is first year when minimum data requirements were met for Control; 2004 (78) for Spawner.
save yr
lb1(1,:) = DataQuadParameters(yr,1,:,2); % Year-specific parameter bounds applied.
ub1(1,:) = DataQuadParameters(yr,1,:,3);
lb2(1,:) = DataQuadParameters(yr,1,:,4);
ub2(1,:) = DataQuadParameters(yr,1,:,5);
Year(yr) "I like to see data printed out from time to time, to make sure that the code is running properly"
for kk = 1:fits
tic

RandVect = rand(1,Paras);
x = lb1 + ((ub1-