How to perform dynamic optimization for a nonlinear discrete optimization problem with nonlinear constraints, using non-linear solvers like SNOPT?

Question

I am new to the field of optimization and I need help in the following optimization problem. I have tried to solve it using normal coding to make sure that I got he correct results. However, the results I got are different and I am not sure my way of analysis is correct or not.

This is a short description of the problem: the objective function shown in the picture is used to find the optimal temperature of the insulating system that minimizes the total cost over a given horizon.

This image provides the mathematical description of the objective function and the constraints

The data of the problems are as follow:

1-

Problem data:

A=1.07×10^8

h=1

T_ref=87.5

N=20

p1=0.001;

p2=0.0037; This is the curve I want to obtain

2- Optimization variable: u_t

3- Model type: The model is a nonlinear cost function with non-linear constraints and it is solved using non-linear solver SNOPT.

4-The meaning of the symbols in the objective and constrained functions

• The optimization is performed over a prediction horizon of N years.
• T_ref is The reference temperature.
• Represent the degree of polymerization in the kth year.
• X_DP Represents the temperature of the insulating system in the kth year.
• h is the time step (1 year) of the discrete-time model.
• R is the ratio of the load loss at the rated load to the no-load loss.
• E is the activation energy.
• A is the pre-exponential constant.
• beta is a linear coefficient representing the cost due to the decrement of the temperature.

I have developed the source code in MATLAB, this code is used to check if my analysis is correct or not.

I have tried to initialize the Ut value in its increasing or decreasing states so that I can have the curves similar to the original one.

This is the curve I obtained

I have tried to simulate the problem using conventional coding without optimization and I got the figure shown above.

close all; clear all;
h=1;
N=20;
a=250;
R=8.314;
A=1.07*10^8;
E=111000;
Tref=87.5;
p1=0.0019; 
p2=0.0037;
p3=0.0037;
Utt=[80,80.7894736842105,81.5789473684211,82.3684210526316,83.1578947368421,... % The value of Utt given here represent the temperature increment over a predictive horizon.
     83.9473684210526,84.7368421052632,85.5263157894737,86.3157894736842,...
     87.1052631578947,87.8947368421053,88.6842105263158,89.4736842105263,...
     90.2631578947369,91.0526315789474,91.8421052631579,92.6315789473684,...
     93.4210526315790,94.2105263157895,95];
 
Utt1 = [95,94.2105263157895,93.4210526315790,92.6315789473684,91.8421052631579,... % The value of Utt1 given here represent the temperature decrement over a predictive horizon.
    91.0526315789474,90.2631578947369,89.4736842105263,88.6842105263158,...
    87.8947368421053,87.1052631578947,86.3157894736842,85.5263157894737,...
    84.7368421052632,83.9473684210526,83.1578947368421,82.3684210526316,...
    81.5789473684211,80.7894736842105,80]; 
Ut1=zeros(1,N);
Ut2=zeros(1,N); 
 
 
Xdp =zeros(N,N);
Xdp(1,1)=1000;

Xdp1 =zeros(N,N);
Xdp1(1,1)=1000;
 
 
   
    for L=1:N-1
     for k=1:N-1
        %vt(k+L)=Ut(k-L+1);
        Xdq(k+1,L) =(1/Xdp(k,L))+A*exp((-1*E)/(R*(Utt(k)+273)))*24*365*h;
        Xdp(k+1,L)=1/(Xdq(k+1,L));  
        Xdp(k,L+1)=1/(Xdq(k+1,L));  
        Xdq1(k+1,L) =(1/Xdp1(k,L))+A*exp((-1*E)/(R*(Utt1(k)+273)))*24*365*h;
        Xdp1(k+1,L)=1/(Xdq1(k+1,L));  
        Xdp1(k,L+1)=1/(Xdq1(k+1,L));  
     end   
        
    end
 % MATLAB code
 
for j =1:N-1
    Ut1(j)= -p1*(Utt(j)-Tref);
    Ut2(j)= -p2*(Utt1(j)-Tref);
        
end
sum00=sum(Ut1);
sum01=sum(Ut2);
 
X1=1./Xdp(:,1);
Xf=1./Xdp(:,20);
Total= table(X1,Xf);
Tdiff =a*(Total.Xf-Total.X1);

X22=1./Xdp1(:,1);
X2f=1./Xdp1(:,20);
Total22= table(X22,X2f);
Tdiff22 =a*(Total22.X2f-Total22.X22);

obj=(sum00+(Tdiff));
ob1 = min(obj);
obj2=sum01+Tdiff22;
ob2 = min(obj2);
 
 plot(Utt,obj,'-o');
 hold on 
 plot(Utt1,obj)