fmincon for bounded optimization problem

Hello,

I am trying to solve a nonlinear optimization problem using fmincon interior point method. Originally my problem formulation does not have bounds on the decision variable, and when i try to run it without the bounds then it takes infintie time and when I run it with bounds then it is much faster. Following are my questions:

1) Why bounds are making the algorithm faster?

2) The final optimal result for the problem is nowhere near the bound, but my lagrange multiplier for the bounds is coming to be non zero, arent they supposed to be zero if the solution is not hitting the bounds?

3) How is the first order optimality criteria defined for interior point method? I saw the documentation but it is not clear to me, is the infinite norm of the grad or some other equation?

I am giving very good initial guess (exact true values) to my problem to make sure it is near the optimal. When I do that, the optimizer is just giving me the intial guess as my final solution which is not possible as I am feeding noisy data to my problem.

Thanks in advance

[John D'Errico]

Whenever I see somone say the solver returns the original point, I want to ask if you are sure your objective is well posed. Is it differentiable? Even if it is theoretically differentiable, are there numerical problems when working in double precision arithmetic? (You would be amazed at what we sometimes do see.) What exitflag is returned? That you see something strange in the return parameters in terms of Lagrange multipliers may be explained by the optimizer crapping out from some problem with the objective, giving up before it even makes an serious effort. So knowing the exitflag is important.

Since we don't see your objective we cannot know what is happening. That means all we can do is offer broad comments.

[Devyani]

function f = objective_CSTR_LAMLE(decision_var,xhat_EKF,Pk_hat_EKF,u_input,Y, Q_MHE_filter,para_est_MHE,C_mat,Qp_est)
global process_para
xkg(1,1:process_para.N_MHE+1) = decision_var(1:process_para.N_MHE+1,1);
xkg(2,1:process_para.N_MHE+1) = decision_var(process_para.N_MHE+2:end,1);
e_x=(xhat_EKF-xkg(:,1));
arrival = e_x'*Pk_hat_EKF*e_x; %computation of arrival cost
f = arrival;
vk(:,1)=zeros(process_para.n_op,1);
wkmhe=zeros(process_para.n_st,process_para.N_MHE);
%specifying Model parameters;
kref_est=para_est_MHE(1);
a_est=para_est_MHE(2)/10^6;
b_est=para_est_MHE(3);
E_by_R_est=para_est_MHE(4)/10^3;
%Loop for MHE
for i = 1:1:process_para.N_MHE
    X0 = xkg(:,i);
    [t,Xt]=ode45(@(t,x)process_dyn(t,x,kref_est,a_est*10^6,b_est,E_by_R_est*10^3,u_input(:,i),[0;0]),...
        [0 process_para.sampling_time], X0);
    Xf = Xt(length(t),:)';
    wkmhe(:,i)=Xf-xkg(:,i+1);
    vk(:,i+1) = Y(:,i)-C_mat*Xf;
    %% Computation of Q
    
        gamma_p(1,1)= -xkg(1,i)*process_para.sampling_time*exp(-1000*E_by_R_est*(1/xkg(2,i) - 1/process_para.Tref));
        gamma_p(1,2)=xkg(1,i)*kref_est*process_para.sampling_time*exp(-1000*E_by_R_est*(1/xkg(2,i) - 1/process_para.Tref))*(1000/xkg(2,i) - 1000/process_para.Tref);
        gamma_p(1,3)=0;
        gamma_p(1,4)=0;
        gamma_p(2,1)=-(xkg(1,i)*process_para.del_H*process_para.sampling_time*exp(-1000*E_by_R_est*(1/xkg(2,i) - 1/process_para.Tref)))/(process_para.Cp*process_para.rho);
        gamma_p(2,2)=(xkg(1,i)*process_para.del_H*kref_est*process_para.sampling_time*exp(-1000*E_by_R_est*(1/xkg(2,i) - 1/process_para.Tref))*(1000/xkg(2,i) - 1000/process_para.Tref))/(process_para.Cp*process_para.rho);
        gamma_p(2,3)=-process_para.sampling_time*((1000000*u_input(5,i)^(b_est + 1)*(xkg(2,i) - u_input(4,i)))/(process_para.Cp*process_para.V*process_para.rho*(u_input(5,i) + (500000*u_input(5,i)^b_est*a_est)/(process_para.Cpc*process_para.rho_c))) - (500000000000*u_input(5,i)^b_est*u_input(5,i)^(b_est + 1)*a_est*((xkg(2,i) - u_input(4,i))))/(process_para.Cp*process_para.Cpc*process_para.V*process_para.rho*process_para.rho_c*(u_input(5,i) + (500000*u_input(5,i)^b_est*a_est)/(process_para.Cpc*process_para.rho_c))^2));
        gamma_p(2,4)=-process_para.sampling_time*((1000000*u_input(5,i)^(b_est + 1)*a_est*log(u_input(5,i))*(xkg(2,i) - u_input(4,i)))/(process_para.Cp*process_para.V*process_para.rho*(u_input(5,i) + (500000*u_input(5,i)^b_est*a_est)/(process_para.Cpc*process_para.rho_c))) - (500000000000*u_input(5,i)^b_est*u_input(5,i)^(b_est + 1)*a_est^2*log(u_input(5,i))*(xkg(2,i) - u_input(4,i)))/(process_para.Cp*process_para.Cpc*process_para.V*process_para.rho*process_para.rho_c*(u_input(5,i) + (500000*u_input(5,i)^b_est*a_est)/(process_para.Cpc*process_para.rho_c))^2));
        
       
        Qdk = gamma_p*Qp_est*gamma_p'+Q_MHE_filter(:,:,i);
    %%
    %computation of objective function for every loop
    f = f + (wkmhe(:,i)'*inv(Qdk)*wkmhe(:,i))+ (vk(:,i+1)'*process_para.inv_R*vk(:,i+1));
end

end

@John D'ErricoThis is the function that is being minimized. CSTR is a two state differential equation that is being solved with ode45, 'f' is the final obejctive function. Equations are differentiable. I have used 'ScaleProblem' as one of the options in fmincon. exitflag is 2 for all time instants

Also, if there are numerical issues in double precision, then what is the general remedy? Also, if I start my code at a different intial guess, then also the results seems to be similar, also the first oder optimality is generally in order of 10^4, which I suppose is high

[Walter Roberson]

We would need all your code to test.

I notice you are using inv() . If you do not have reason to suspect that Qdk might be singular, then use

wkmhe(:,i)'/Qdk*wkmhe(:,i)

If you do have reason to suspect that Qdk might be singular, then use pinv() instead of inv()

[Devyani]

Okay I am attaching all the files for your reference.

Thanks.

[Devyani]

@Walter Roberson Hello, I have given my code. Also i checked with inv, the Qdk is not singular.

Thanks in advance

[Walter Roberson]

I had to uninstall MATLAB temporarily because of operating system limitations (not related to MATLAB itself.) It may take me a bit of time to recover.

[Devyani]

Okay thank you @Walter Roberson

Answer 1

[Walter Roberson]

Generally speaking, functions without bounds can take indefinite time to minimize if the function has an asymptope

        |    |
    ___/ --v-+
         

where the v marks the minimum. But if the function happens to land on the shoulder to the left then the local gradient slopes away from the center and the minimizer can take indefinite time exploring that left slope.

fmincon is a local optimizer: there is no way for it to know that it should spend time climbing the hill to the center to get to a better opportunity. And not every such hill happens to lead to a global minimum. Furthermore, the global minimum can be an indefinitely small part of the graph -- imagine putting something heavy on a section of stiff rubber that has a very elastic center, then unless you were quite close to the indentation you would get no information that it existed.

[Devyani]

Thank you for the response, that explains why I need bounds. But why the langrange multiplliers are non-zero when the final optimal values are not near the bounds.?

[Walter Roberson]

Sorry, I am not familiar with the theory about the lagrange multipliers. (I read it once, but did not retain it in memory.)

[Devyani]

Okay thank you anyways for solving the first part.:)