Optimization of simulation program result w.r.t. program parameters

[Rubén Roa-Ureta]

Hi:
I need to optimize the final output (a real number) of a long program 
w.r.t. three parameters that determine how the program works. I have 
prepared the smallest toy example I could think of below. My purpose is 
to compare the program's prediction of the final number of individuals 
in a cohort (as a result of the parameters that determine egg 
development and egg mortality in the early life) with the same number 
obtained from another sampling-model process, to make these two numbers 
as close as possible.
I have found information to optimize a function (a mathematical 
expression), for example with samin and bfgsmin, but I am not sure if 
this is the right strategy to optimize the output of a long program. In 
my program, there are a few places where random numbers are called and 
used to create variety among the individuals.
The toy example is the minimum implementation of an individual-based 
model; in this case I have 1000 super-individuals and each one starts 
containing the same starting number of individuals but they differ 
because they have a different potential of temperature-dependent growth 
and wind-dependent mortality.
Shall I treat this problem as a function optimization problem, turning 
my program into a function to be called by the optimizer?
Toy example follows. Real program takes about 3 minutes to finish on a 
Windows Vista Core Duo at 1.73 Ghz and 1 GB RAM.
Thanks in advance
Ruben

%par.par, text file with five real numbers
par = load('par.par');
N0 = par(1); % start number of eggs, very large, ca. 10^12
b = par(2); % development parameter 1
c = par(3); % development parameter 2
d = par(4); % egg mortality parameter
e = par(5); % juvenile mortality parameter
f = par(6); % day when egg turn to juvenile
g = par(7); % final day of simulation
Nt = par(8); % number of individuals in final day, known from another 
model-sampling process
%data, text file with two numeric columns and n rows, one column for 
temperature, one for wind speed.
data = load('data.dat');
temp = data(:,1);
wind = data(:,2);
ind = 100000; number of super-individuals
winter = zeros(ind,6); %matrix containing the information about 
super-individuals, overwritten at each time step.
%Column 1: Start number;
%Column 2: Temp-dependent development;
%Column 3: Surviving numbers;
for i=1:ind
  winter(i,1) = a/ind; % starting number of individuals. Does not change 
with time
  winter(i,2) = random growth potential, uniform number between 0.5 and 1.5;
  winter(i,3) = random susceptibility to wind-related mortality, normal 
between 0.75 and 1.25;
  winter(i,4) = development(a,b)*winter(i,2); turns to 1 when 
development of egg into juvenile is completed
  winter(i,5) = day when egg turns into juvenile; when winter(i,4) turns 
to 1
  winter(i,6) = number of individuals, at first increase due to eggs 
transforming into juveniles (but controlled by wind-dependent egg       
mortality), then decrease due to egg and juvenile tomortality
end
for day=1:g
  winter(i,4) = degree of development dependent on daily temperature an 
individual growth potential
  winter(i,5) = written with index 'day'  winter (i,4) turns to 1
  winter(i,6) = gets updated with temperature-dependent transformation 
of eggs into juveniles and wind-
  dependent mortality of eggs, until day f, and later, with juvenile 
mortality rate (independent of wind)
end
compute number of individuals in g, final day, Nt-program.
compute squared difference between Nt-program and Nt observed from 
another model-sampling process
optimize squared difference w.r.t. parameters b,c,d