Optimization in a for loop

Hi,

I have the following code which is a simiplification of a more complex one. The code contains an optimization problem (prob) that finds the optimum values of U1 and U2 to minimize the Root-Mean-Square-Error (prob.Objective) between the simulated (Heat_flow) and measured (M) data. As you can see, the code works and actually finds the optimum solutions. However, those solutions vary for each iteration of the loop.

What I would like to do, is to have only two optimum solutions (for U1 and U2, respectively) that do not vary at each iteration. I tried to introduce a constraint (prob.Constraints.U1=eq(U1,U1(1))), but the solver finds a solution for the first iteration and then applies it to the following ones.

How can I do it?

Thank you!

%Input data (outdoor temperature and areas)
T_out=[18 16 14 15 19 14 17];
Area_1=16;
Area_2=19;

%Measured data
M=[79 84 64 76 86 56 76]';


%Definition of the optimization problem
prob=optimproblem("Description","Esempio");

%Definition of the optimization variables (U-values)
U1 = optimvar('U1',7,1,"LowerBound",0);
U2 = optimvar('U2',7,1,"LowerBound",0);

%Calculation of the heat transfer coefficient
Heat_transfer=U1*Area_1+U2*Area_2;

%Loop for perfoming the calculation for each measurement
Heat_flow=[];

for i=1:max(size(T_out))

    Heat_flow=Heat_transfer*(20-T_out(i));

end

%Definition of the objective function (minimization of RMSE) 
prob.Objective=((sum((M-Heat_flow).^2))/7)^(1/2);

%Definition of the initial guess (not a linear problem)
initialGuess.U1 = [1 1 1 1 1 1 1];
initialGuess.U2 = [1 1 1 1 1 1 1]; 


%Solution of the problem
[sol,opt]=solve(prob,initialGuess)

[Sam Chak]

Hi @Andrea Costantino

The optimization problem is not formulated correctly, specifically the objective function. Can you describe the objective function, both looped and without loop?

% Input data (outdoor temperature and areas)
T_out  = [18 16 14 15 19 14 17];
Area_1 = 16;
Area_2 = 19;

% Measured data
M      = [79 84 64 76 86 56 76]';

% Definition of the optimization problem
prob   = optimproblem("Description", "Esempio");

% Definition of the optimization variables (U-values)
U1     = optimvar('U1', 7, 1, "LowerBound", 0);
U2     = optimvar('U2', 7, 1, "LowerBound", 0);

% Calculation of the heat transfer coefficient
Heat_transfer = U1*Area_1 + U2*Area_2;

% Loop for perfoming the calculation for each measurement
Heat_flow = [];

% Disable this the loop
% for i = 1:max(size(T_out))  % loop till i = 7
%     Heat_flow = Heat_transfer*(20 - T_out(i));
% end

% to test this Heat flow:
Heat_flow = Heat_transfer*(20 - T_out(7));

% Definition of the objective function (minimization of RMSE) 
prob.Objective = (sum((M - Heat_flow).^2)/7)^(1/2);

show(prob)

  OptimizationProblem : Esempio

	Solve for:
       U1, U2

	minimize :
       (sum((extraParams{1} - (((U1 .* 16) + (U2 .* 19)) .* 3)).^2) ./ 7).^0.5

         extraParams{1}:

         79
         84
         64
         76
         86
         56
         76





	variable bounds:
       0 <= U1(1)
       0 <= U1(2)
       0 <= U1(3)
       0 <= U1(4)
       0 <= U1(5)
       0 <= U1(6)
       0 <= U1(7)

       0 <= U2(1)
       0 <= U2(2)
       0 <= U2(3)
       0 <= U2(4)
       0 <= U2(5)
       0 <= U2(6)
       0 <= U2(7)

%Definition of the initial guess (not a linear problem)
initialGuess.U1 = [1 1 1 1 1 1 1];
initialGuess.U2 = [1 1 1 1 1 1 1]; 

%Solution of the problem
[sol, opt] = solve(prob, initialGuess)

Solving problem using fmincon.

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

sol = struct with fields:
    U1: [7×1 double]
    U2: [7×1 double]

opt = 6.3137e-09

[Matt J]

However, those solutions vary for each iteration of the loop.

That's not apparent to me. The code you have presented does not show multiple optimization problems in a loop. In your code, solve() is called only once.

[Matt J]

It would be advisable to reformulate your objective in purely quadratic form. Then, more specialized solvers like lsqlin can be used.

prob.Objective= sum( (M-Heat_flow).^2 ); %equivalent to what you had before.

[Andrea Costantino]

Yes, solve() is called once, but the optimization variables are 7x1. I had to define them as 7x1 because 7 iterations are present. To create Heat_flow in the for loop, a new optimization expression has to be created in each iteration since it is written as a symbolic expression.

Answer 1

[Matt J]

Perhaps this is what you wanted. If so, it is a rather ill-posed problem, because your C matrix is only rank 1.

%Input data (outdoor temperature and areas)
T_out=[18 16 14 15 19 14 17];
Area_1=16;
Area_2=19;

%Measured data
M=[79 84 64 76 86 56 76]';


%Definition of the optimization problem
prob=optimproblem("Description","Esempio");

%Definition of the optimization variables (U-values)
C=(20-T_out(:))*[Area_1,Area_2];
U = optimvar('U',2,1,"LowerBound",0);


%Definition of the objective function (minimization of RMSE) 
prob.Objective=sum(C*U-M).^2;

%Solution of the problem
[sol,opt]=solve(prob)

Solving problem using lsqlin.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = struct with fields:
    U: [2×1 double]

opt = 1.1953e-22

[Andrea Costantino]

Thank you Matt J!

I understand the process. Basically, you skip the loop by using vectors operation. Unfortunately, in my real code I cannot do it because each iteration represents a step and the elements of the vector are computated from the previous iteration.

[Matt J]

That's irrelevant. Each row C(i,:) represents the coefficients for a particular step. I could have used a loop to create C in a way that depends in some way on previous steps, e.g.,

C=nan( numel(T_out), 2);

for i=1:height(C)
    
    C(i,:)=(20-T_out(i))([Area1,Area2]);
    
    if i>1
        C(i,:)=C(i,:) + C(i-1,:);
    end
    
end