Optimization problem

Specific domains Optimization
1

Dear,

I always struggle to understand code for optimization problems.

Hence I kindly ask for your help with following task:

The task is to optimize a set of calibration coefficients, such that a calculated value is as close as

possible to a reference…in other words: The distance between the calculated value and the reference value is minimal.

There are some start values for the calibration coefficients and and also some upper and lower bounds for them.

Also a stop criteria is defined and we want to introduce a max number of iteration steps

// start values for calib coeff
calibCoeff_0 = [ t1_0, t2_0, p1_0, p2_0, p3_0];

// upper bounds for calib coeff
UBDS = [ t1_up, t2_up, p1_up, p2_up, p3_up,];

// lower bounds for calib coeff
LBDS = [ t1_lo, t2_lo, p1_lo, p2_lo, p3_lo,];

// stop criteria ... end optimization if the result of myfunc() < STOP_CRIT
STOP_CRIT = 1*10^(-6)

// maximal iteration steps
MAX_ITER_STEP = 6000;

// library function to calculate P & T values based on calibCoeff and raw values

[P, T] = calcValues(calibCoeff, p_Raw, t_Raw);

I think a function is necessary that compares the result using the actual coefficients against the reference value

// function to compare result with actual calibCoeff against reference value
// "res" shall be minimized

function [f, g, ind] = myfunc([pRaw tRaw refVal], ind, parameters)

   // assumption [pRaw tRaw refVal] = x from help page ...is this correct(?)
                

    // calculate P and T based on current calibCoeff (parameters)
    [P, T] = calcValues(parameters, pRaw, tRaw); 
    
    // calcluate difference from reference
    res = (refVal - P)^2

    // assumptions ... is that correct(?)
    f = P;
    g = res;


endfunction

Now: I GUESS that myfunc() is what is normally called the cost function. As this needs a set of parameters I would follow the approach from the help page of using a list

// Store the parameters
parameters = tlist ([
"T_MYPARAMS"
"t1"
"t2"
"p1"
"p2"
"p3"
]);

parameters.t1 = t1;
parameters.t2 = t2;
parameters.p1 = p1;
parameters.p2 = p2;
parameters.p3 = p3;

costf = list (myfunc, parameters);

I GUESS the optim call would look similar to this….

[fopt, xopt] = optim (costf, calibCoeff, "b",  LBDS, UBDS, calibCoeff_0, MAX_ITER_STEP, STOP_CRIT)

However:

As you may have realized I am a bit confused and any help would be appreciated. e.g.: would xopt be the set of optimized calibCoeff (parameters) ?

Thank you,

Philipp

2

I think you messed up with parameters (pRawn, tRaw, refVal) vs. unknowns (calibration coefficients) of your optimization problem. The cost function should start like this (I supposed that P is actually a vector, hence the calculation of the error should sum the squared residues)

function [f, g, ind] = myfunc(calibCoeff, ind, pRaw, tRaw, refVal)
    [P, T] = calcValues(calibCoeff, pRaw, tRaw); 
    f = sum((refVal - P).^2);
    ...
end

and costf should be defined as (you can replace the parameters sequence by a struct or a tlist, but for 3 parameters it is simpler to pass them directly)

costf = list(myFunc, pRaw, tRaw, refVal);

However (see the continuation dots in myfunc above), you are in trouble because you have to provide the gradient of myFunc with respect to califCoeff. If you master what goes on in calcValues you can obtain the gradient by hand, but if this is not the case you have to approximate it by using NDcost (see NDcost - Generic external for optim computing gradient using finite differences), this can be done by redefining costf like this

costf = list(NDCost, myFunc, pRaw, tRaw, refVal);

The optim call will give you the optimized values of calibration coefficients:

[fopt, calibCoeff] = optim (costf, "b",  LBDS, UBDS, calibCoeff_0, "ar", MAX_FCALL, MAX_ITER_STEP, STOP_CRIT)

where MAX_FCALL is the maximum number of cost function call.

If using toolboxes is OK for you, fmincon (see ATOMS : Fmincon details) uses a more robust solver (interior-point method) and has a better user, logging and tuning interface.

S.

3

Hello Stéphane,

Thanks for your input.

To clarify:

calibCoeff a set of sensor calibration coefficients [t1, t2, p1, p2, p3] initially provided by manufacturer

For each sensor one could write:

calibCoeff.t1 = t1    // above I used parameter.t1 to be in line
calibCoeff.t2 = t2    // with the help page
calibCoeff.p1 = p1
calibCoeff.p2 = p2
calibCoeff.p3 = p3

calcValues()is a self written function to calculate P-T-Values (Pressure & Temp) from sensor raw values (pRaw, tRaw) by use of calibration coefficients [t1, t2, p1, p2, p3]

calcValues() works with scalars and / or vectors

gradient calculation Since f can be calculated , would it be enough to track the difference between f(i) and f(i-1) such as:

gradient = f(i) - f(i-1) // i = a single step within the optimization process

Then: if gradient < STOP_CRIT –> stop optimization

Again: My understanding is most likely insufficient.

fmincon

using toolboxes would be OK … can fmincon be updated for Scilab 2025.1.0 or higher?

4

You can give a first try without providing the gradient. The actual algorithm computing the true gradient is completely dependent of the code of your calcValues function so without seing the formulas (this is easier than working on the code implementing the formulas).

However, a thing that I missed is that your problem is a least squares problem. You may use leastsq(see leastsq - Solves non-linear least squares problems), which can work without providing the gradient (but relies internally on optimand NDCost). In that case the user function should provide the residual vector (and not its squared norm). So the code will be more compact:

function r = myresidual(calibCoeff, pRaw, tRaw, refVal)
    [P, T] = calcValues(calibCoeff, pRaw, tRaw); 
    r = refVal - P;
end
resid = list(myResidual, pRaw, tRaw, refVal);
[fopt, calibCoeff] = leastsq (resid, "b",  LBDS, UBDS, calibCoeff_0, "ar", MAX_FCALL, MAX_ITER_STEP, STOP_CRIT)
5

Thank you very much.

In general this seems to work just fine.

But sometimes I discover is that the optimizer output results that does not fit.

E.g.: Using the optimized coefficients the the calculated pressures will be way off of reference.

In this cases I need to play with LBDS and UBDS … that seems to help.