Solving an optimization problem in fmincon with nonlinear constraints

Hi All,

I have the following dynamical systems.

Equation 1:

Equation 2:

Equation 1 represents the exact dynamics of a system and equation 2 is the approximate dynamics that will give the same time course profiles as equation 1, upon optimization. For the sake of understanding, let us assume eq(1) gives experimental values and equation 2 are the predicted values.

The objective function defined has a cost function that minimizes the difference between state variables ϕ and , by optimizing parameter which are the control variables.

I'm trying to solve this as a parameter estimation problem with non-linear equality constraints .The non-linear constraints have been written by discretizing equation (2) using trapezoidal collocation.

I have defined function [c ceq] = defects( ) that accepts as input argument and the and d values necessary for setting up the constraints are obtained by calling a nested function defined within defects` .This nested function accepts input arguments ode15s(@(t,) fun(t,, )`.

I've implemented the above , a few iterations are carried out by fmincon and there is a dimension error that pops up. I am not sure why this occurs after a few iterations and not at the first iteration itself.

Error using cellfun
All of the input arguments must be of the same size and shape.
Previous inputs had size 1 in dimension 1. Input #4 has size 51

Error in cse_03_19_20/defects (line 117)
f = cellfun(@predicted,num2cell(t),num2cell(phi_tilde,2),num2cell(Dhat_mat,2),'UniformOutput',false);

Error in barrier

Error in fmincon (line 824)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in cse_03_19_20 (line 9)
Dhat= fmincon(@objfun,Dhat0,[],[],[],[],[],[],@defects)

Could someone help me in finding out what leads to dimension error after a few iterations? Please find the code attached

EDIT: I've updated the file attached and also used dbstop if error as suggestd below

Now the error that occurs is

Matrix dimensions must agree.

Error in cse_03_19_20/objfun (line 31)
    f = (phi - phi_tilde).*(phi - phi_tilde);

Error in barrier

Error in fmincon (line 824)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in cse_03_19_20 (line 10)
Dhat = fmincon(@objfun,Dhat0,[],[],[],[],[],[],[], opts)

But the same code gives optimal solution when lsqnonlin is used instead of fmincon.

Any suggestions on why fmincon fails to give optimal solution will be helpful

[Deepa Maheshvare]

Any suggestions?

[Ameer Hamza]

I tried running your code. I didn't get any such error, but the fmicon was also not able to reach an optimal solution. To track down the cause of this error, you can type this in the command window

dbstop if error

and then run your code. This will tell MATLAB to halt execution whenever the error occurs and show you the line of code. You can then check the values of the variable and try to run the line manually to see what is causing the error.

[Deepa Maheshvare]

Could you please check my edit?

[Ameer Hamza]

I tried the code with both fmincon and lsqnonlin without any error. Both give comparable solutions. Have you attached the exactly same code you are running?

fmincon:

Dhat =

   1.0e+03 *

    4.9760
    4.8535
    5.3041
    5.0159
    5.1509
    4.3784
    4.7378
    7.9296
    5.0528

lsqnonlin:

Dhat =

   1.0e+03 *

    5.0608
    4.9864
    5.0414
    5.1345
    5.1632
    5.0217
    4.5965
    3.7518
    2.2228

[Deepa Maheshvare]

Thanks for the response. I am not really sure what's happening

Here are my results,

lsqnonlin

Dhat =

   1.0e+03 *

    5.0005
    5.0002
    5.0002
    5.0002
    4.9997
    4.9969
    4.9811
    4.8900
    4.9051

fmincon % simulationg without nonlinear constraints defined in my original post

Matrix dimensions must agree.

Error in cse_03_19_20/objfun (line 37)
    f = (phi - phi_tilde).*(phi - phi_tilde);

Error in barrier

Error in fmincon (line 824)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in cse_03_19_20 (line 14)
   Dhat = fmincon(@objfun,Dhat0,[],[],[],[],[],[],[], opts)

chnage to make in attached file: line5

solver = struct('lsq',false); while running fmincon and tru while running lsqnonlin

[Deepa Maheshvare]

Please check here too , I posted this as a separate question since I'm testing the problem with no constraints now

[Deepa Maheshvare]

Hi Ameer,

Did you have a chance to look at the updated file? I hope you could oberve the same results that I reported above

[Ameer Hamza]

@Deepa, i again got the same result with fmincom and lsqnonlin. However, Walter also reported the same error as you. I wonder if this related to the version of MATLAB. I tried it on R2020a.

[Walter Roberson]

I agree, I do not see the problem in R2020a, but do see it in R2019b.

[Deepa Maheshvare]

Thanks Walter. I will try to install R2020a and check

[Deepa Maheshvare]

@Ameer and @Walter. MATLAB version 2020a is not available at my institute yet. Please let me know if this works in lower version of MATLAB (other than 2019b that I am currently using.

Or, is there other alternative?

[Ameer Hamza]

Deepa, as Walter mentioned that it works fine for R2020a but fails for R2019b. This shows that Mathworks has improved the fmincon function with their new release. I am not sure how to solve the problem in R2019b, but as Walter mentioned, when Dhat becomes large, the gradient estimation fails. So you somehow need to reformulate the problem to avoid the infinities.

[Walter Roberson]

Same problem in R2018b.

[Deepa Maheshvare]

Thanks for the update. Can I get a trial version of 2020a?

It is not really straightforward to reformulate the problem :( . Dhat goes into differential equations.

[Walter Roberson]

https://www.mathworks.com/campaigns/products/trials.html

[Deepa Maheshvare]

Thanks Walter!

From what is mentioned here https://in.mathworks.com/products/matlab-online.html#license-types , I cannot use MATLAB online with trial license?

[Deepa Maheshvare]

@Ameer and @Walter Using 2020a resolved the issue. However, I would like to know why there is a large difference in the optimal values reported by lsqnonlin and fmincon

[Ameer Hamza]

fmincon and lsqnonlin use diifferent optimization algorithm. You can calculate the value of objective function at both solutions to see if there is a significant difference.

[Deepa Maheshvare]

@Ameer

Thank you. I checked the default algorithm used by both fmincon and lsqnonlin

fmincon

algorithm : interior point (default)

lsqnonlin:

algorithm : trust region reflective(default)

I tried to use trust region reflective in fmincon (because trust region reflective used in lsqnonlin gives more accurate predictions). But I couldn't. In the documentation on `Choosing the algorithm`

it is mentioned "'trust-region-reflective' requires you to provide a gradient"

I am not sure how to provide gradient for my objective function.

For instance, no I have 10 variables. Should I compute the partial derivative of the objective function wrt 10 variables?

If yes, this will be a vector ]

Is there a way to do this in MATLAB automatically? And how should these gradients be evaluated?

Any explanation in this direction will be really helpful.

[Ameer Hamza]

Theoretically, you can write down your objective function and analytically find its gradient by using a little bit of knowledge about multivariate calculus linear algebra. However, I don't think it will significantly improve the performance since fmincon and lsqnonlin are already given the solution of the same order of magnitude.

Answer 1

[Walter Roberson]

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In ode15s (line 589)

In cse_03_19_20/objfun (line 32)

In barrier

In fmincon (line 824)

In cse_03_19_20 (line 14)

(a whole bunch of those, then)

Warning: Failure at t=2.205524e-02.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(5.551115e-17) at time t. 
> In ode15s (line 668)
  In cse_03_19_20/objfun (line 32)
  In barrier
  In fmincon (line 824)
  In cse_03_19_20 (line 14) 
Matrix dimensions must agree.

Error in cse_03_19_20/objfun (line 37)
    f = (phi - phi_tilde).*(phi - phi_tilde);

Error in barrier

Error in fmincon (line 824)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in cse_03_19_20 (line 14)
   Dhat = fmincon(@objfun,Dhat0,[],[],[],[],[],[],[], opts)
 
37      f = (phi - phi_tilde).*(phi - phi_tilde);

Notice the failure at t=2.205524e-02. Because of that, the ode15s call is returning early, before it has processed all 51 timestamps that you requested output at. The first call for exact solutions did all 51, but the second call for approximate only gets 3 of them. When you do phi - phi_tilde under that circumstance, you are subtracting arrays of different sizes.

[Walter Roberson]

When Dhat gets large enough such as

[5409.35819159269 4457.86895613171 4338.14787977186 5841.61446846487 ...
 5286.48060962025 265.616463594906 -4898.82021376954 -11384.6129334742 ...
 -7092.42118240817] .'

then the integration grows without bound, such as up to

[5 1.01274196457643e+297 9.61416483594012e+297 ...
 9.03733004653503e+298 6.52402843345009e+299 5.2783698527968e+300 ...
 7.42243004636899e+302 -4.21470318859266e+303 5.66646852034728e+303 ...
 -4.40008055621133e+303] .'

At that point, the A*phi_hat overflows to Inf on one of the entries. Once that happens, you get Inf propagating, and that quickly leads to infinities minus infinities, which gives you NaN entries. The gradient estimation then fails, leading to ode15s exiting early.

[Ameer Hamza]

That's really strange. I ran it several times with fmincon, and it always reaches the same solution.

Local 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.

<stopping criteria details>

Dhat =

   1.0e+03 *

    4.9760
    4.8535
    5.3041
    5.0159
    5.1509
    4.3784
    4.7378
    7.9296
    5.0528

Can you confirm that you tried it on R2020a or some earlier release?

[Deepa Maheshvare]

Hi Ameer,

I am trying it on R2019b.

[Walter Roberson]

Notice that the vector I indicated as giving a problem, Dhat =

[5409.35819159269 4457.86895613171 4338.14787977186 5841.61446846487 ...
 5286.48060962025 265.616463594906 -4898.82021376954 -11384.6129334742 ...
 -7092.42118240817] .'

has negative entries, but the R2020a result only has positive entries. A speculation would be that negative entries beyond some threshold lead to a feedback loop that gives overflow.

Is there justification for constraining the DHat values?

[Deepa Maheshvare]

Sorry for the delay in response.

Equations presented in my original post:

Equation 1:

Equation 2:

In the code I am basically solving the same system in equation 1 and 2. In eqn 1 , I have used D to simulate the ode . In equation 2, I'm changing D to D/10 (which I represent as . Now , I am trying to see if the optimization solver can generate back D by reducing the difference in time course data generated by solving equation 1 and 2. Using lsqnonlin and fmincon , I have tried to solve this as an unconstrained optimization problem .

[Walter Roberson]

I have a custom minimizer working on the problem. It has been running several days so-far and has not completed even one pass yet. I would interrupt it and tell it to use a less refined search, but then I would lose those days of work...

[Deepa Maheshvare]

The values obtained from lsqnonlin and fmincon are quite different although they lie in the same order of magnitude.

fmincon:

Dhat =

          4964.15968380843
          4844.55968898685
          5338.31296991717
          5033.10965791698
          5141.84022911585
          4323.44217676052
           4791.4163677027
          8052.78665412938
          5128.19019018659


output = 

  struct with fields:

         iterations: 116
          funcCount: 1174
    constrviolation: 0
           stepsize: 3.95654181075794
          algorithm: 'interior-point'
      firstorderopt: 9.96259102784739e-07

lsqnonlin:

Dhat =

          5000.47884987158
          5000.21574114392
          5000.22782217235
          5000.17680445071
          4999.71112597238
          4996.88552944583
          4981.07010047556
          4890.07762209202
           4905.6437322249


output = 

  struct with fields:

    firstorderopt: 6.02617504535526e-07
       iterations: 13
        funcCount: 140
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 885.613142141269

lsqnonlin is giving more accurate results in less number of iterations. The output close to the values set for D , in equation 1. I tried to use the same algorithm (i.e trust-region-reflective) in `fmincon` but it requires gradients specified by the user. At this point, I am not sure how that has to be implemented. So, I haven't tried this yet.

[Walter Roberson]

After numerous days of my computer working on the problem, it turns out that there is a perfect fit (exact zero return) when Dhat is all exactly 5000.

[Deepa Maheshvare]

Hi Walter,

Thanks a lot for letting me know. That's a great news! Both lsqnonlin and interior-point algorithm did not give all 5000. I'm kind of stuck here.

Could you please let me know the algorithm that has been used in your custom solver? I'd happy to know if the code is already available in any of the public repositories

Was there any initialization strategy that has been used? Is there an in-build way to compute the soluton of odes inside the objective function ? Lately, I tried to solve the same problem in python-GEKKO and I am facing the same problem in GEKKO too.

[Walter Roberson]

My custom algorithm uses a mix of strategies that in part involves a lot of brute force searches to get into the right area. That got me a number of positions very close to all 5000, so afterwards I had it calculate all 5000 and sure enough got a perfect fit.

[Deepa Maheshvare]

Thanks a lot for sharing details. Could you please offer some suggestions on how to proceed further?

[Walter Roberson]

There isn't any "further" to proceed to. You pass all 5000's into your function and you see that you get exactly 0 back. You examine your function and see that you are effectively calculating a norm, and you observe that norms can never be negative, so you realize that your function can never return negative and so 0 is the optimal for it.

If your goal is to optimize this function, then you take the 5000's and you are done with it.

If your goal is to have this function as a sample problem to demonstrate how to get the optimization tools to find global minima, then you also need to observe that the input options on tolerances make a notable difference into when to give up. If you did not know that a perfect fit was possible, and you were getting results from your function that were on the order of 1e-21 with all the values near to 5000, then how would you know whether to give up or not? If you want the numeric optimizers to prove that they have reached a global minima, then you need to give up on that idea -- symbolic calculations are the only hope to prove that you have reached a global minima.