optimization of temperature profile

Hello.

I wish to use the optimization solver to solve an optimization problem. i have a vector of Temperatures (T). with the initial condition

T(1) = T_min = -30

in the code it's cons 3

i also have constrain on the gradeint of the temperature, in deg/sec:

0.1 <= (T(i+1)-T(i))/dt <= 5

in the code it's cons 1 and 2

i wish to minimize the following cost function

sum( (T-Ts))^2

where Ts is some temperature i would like to be closer to them ( in the code it's -30, -5 35 and 60)

i also would like it to be eqully spred, in the code it's cons4.

i tried play with the options but i can't seem to figure it out beacuse the result it not finding an optimom that aplly to all of the constrains (even the initial temperature is not close to -30 for example), any ideas?

here my code:

T_max      = 60;     % C
T_min      = -30;    % C
T_dot_max  = 5/60;   % C/sec
T_dot_min  = 0.1/60; % C/sec
Time       = 10000;   % sec
dt         = 10;
N          = Time/dt;
bins       = 20;

Ts = [-30 -5 35 60];
f  = [1/4 1/4 1/4 1/4];
%% Optimization problem
prob  = optimproblem;
T     = optimvar('T',N,'LowerBound',T_min,'UpperBound',T_max);
t     = 0:dt:Time;
%setting initial point - change here
x0.T     = T_min ;

%% Constrains
for i=1:N-1
  cons1(i) = (T(i+1)-T(i))/dt >= T_dot_min;  
  cons2(i) = (T(i+1)-T(i))/dt <= T_dot_max;
end
cons3 = T(1) == T_min;


for j =1:length(Ts)
    for i=1:N/4
        SF((j-1)*(N/4) + i) = ((T(i)-Ts(j))^2)*dt^2;
    end
end

for j =1:length(Ts)
    SFF(j) = sum(SF((j-1)*N/length(Ts) + 1:(j)*N/length(Ts)));
end


for j =1:length(Ts)-1
    cons4(j) =  SFF(j) == SFF(j+1);
end
prob.Constraints.cons1 = cons1;
prob.Constraints.cons2 = cons2;
prob.Constraints.cons3 = cons3;
prob.Constraints.cons4 = cons4;
%% Cost function
%cost func

one = ones(N,1);
for j =1:length(Ts)
    for i=1:N/4
        SUM((j-1)*(N/4) + i) = ((T(i)-Ts(j))^2)*dt^2;
    end
end

prob.Objective  =  sum(SUM(:));

for j =1:length(Ts)
    for i=1:N/4
        x0.T((j-1)*(N/4) + i) = Ts(j);
    end
end
x0.T = T_min + (T_max-T_min)*(1:N)./N;
options = optimoptions('fmincon','MaxIterations',8e10,'Algorithm','sqp');
options.MaxFunctionEvaluations   = 6e+20;
options.ConstraintTolerance      = 1e-20;
options.StepTolerance            = 1e-20;
options.FunctionTolerance        = 1e-20;
options.OptimalityTolerance      = 1e-20;
options.StepTolerance             = 1e-20;
[sol,fval,exitflag,output] = solve(prob,x0,'Solver','fmincon','Options',options);
disp(sol);