how to create constraints in fmincon

Hi all,

I have a dataset with two variables x and y as follows:

10 0.006 11 0.017 12 0.026 ...

50 0.003 51 0.002 52 0.001

I am trying to fit it with a nonlinear function with 4 parameters a, b, c, and d, of which c and d should be 10<c<52 and 10<d<52. The estimated y should be greater than 0. How can I create such constraints in fmincon? I was able to get the estimate using fminsearch. But many time, c and d and y do not meet the criterion.

lb and ub options in fminsearch are for y, not for parameters (right?) I am new with Matlab.

Thank you for warm support!

Rdu

[Yu Jiang]

• In the function fmincon (see documentation), the input arguments lb and ub are in fact the lower and upper bound for the decision variables a, b, c, d, but not for the output y. If you need to keep y as a non-negative number, you may need to impose that constraint by choosing appropriate coefficient matrices Aeq and Beq.
• It would be helpful if you can provide more details regarding the nonlinear function you are using and the objective function that you are minimizing by using fmincon.
• I also wonder if you have tried to use the Curve Fitting Toolbox (see documentation) , which seems to be more suitable to be used here than fmincon.

-Yu

Answer 1

[Yu Jiang]

This should be an optimization problem with nonlinear constraints, and the problem should be able to be solved using fmincon (See Documentation).

Before using fmincon, two functions need to be defined. One is the objective function, one is the nonlinear constraint function.

Let us denote the nonlinear function as foo(x,a,b,c,d), which is the mixture Beta distribution as you mentioned previously. Then, the objective function can be defined as

% Objective function
function err = myObj(p,x,y)
 a = p(1);
 b = p(2);
 c = p(3);
 d = p(4);
 err = 0;
  for i = 1:length(x)
    err = (foo(x(i) , a, b, c, d) - y(i))^2;
  end
end
% End the of the objective function

Save the above two functions as two separate MATLAB files of which the filenames are the same as the function names.

Now, the nonlinear constraint also needs to be defined. Since it is required that f(x,a,b,c,d)>=0 for all x, a way to defined the nonlinear constraint is

% Constraint function
function [nl,nl2] = nlcon(p,x)
 a = p(1);
 b = p(2);
 c = p(3);
 d = p(4);
 nl = [];
 for i = 1:length(x)
    nl = [nl; -foo(x(i), a,b,c,d)]; 
 end
 nl2= [];
end
% End of the constraint function

Notice that I put a negative sign before foo. This is because fmincon requires the nonlinear constraint function to be smaller or equal to zero.

Now, in MATLAB command line or a script file, you can try the following steps

Step 1) Give values to x and y, such as

>> x =[10:52]; y = [0.006, 0.017, 0.026, ,..., 0.003, 0.002, 0.001]

Step 2) Specify the upper and lower bounds

>> lb = [10; 10; -inf; -inf]; % Assuming 10 is the lower bound for |a| and |b|
>> ub = [52; 52; inf; inf]; % Assuming 52 is the upper bound for |a| and |b|

Step 3) Specify the initial condition

>> p0 = [5;5;5;5]; % Give some initial guess for a,b,c,d

Step 4) Solve the optimization problem using fmincon

>> xx = fmincon(@(p)myObj(p,x,y),p0,[],[],[],[],lb,ub,@(p)nlcon(p,x))

Note that in the above code I use anonymous functions (See Documentation) to pass x and y. This is because fmincon requires the objective function and the constraint function to only contain the decision variable p as their arguments, while the functions I defined have (p,x,y) or (p,x) as the arguments. To ensure the functions satisfy the requirements of fmincon, I defined to anonymous functions which only depend on p.

-Yu

[Rdu]

Hi Yu,

Many thanks for your generous efforts and help. Your code is easy to understand. However, the code did not let me pass through, mainly from the line of "err=" in myObj. My x and y are listed below.

if true
  % 
x=[12:52];
y=[0.02  0.04    0.06    0.07   0.16  0.25  0.32  0.38  0.43...  
 0.47   0.50  0.51     0.51   0.51  0.50  0.49  0.47  0.45...
 0.43   0.41  0.39    0.36  0.34  0.32  0.29  0.27  0.24...
 0.20   0.17  0.14    0.11  0.09  0.07  0.05  0.04  0.03...
 0.02   0.02   0.01    0.01    0.01];

end

I could not figure out what's wrong with the code. I tried many different ways, all failed. My foo function is as follows

if true
  % function f =foo(x,p)

 a = p(1);
 b = p(2);
 c = p(3);
 d = p(4);

f=(gamma(a+b)./(gamma(a).*gamma(b))).*((d-c).^(-(a+b-1)))... .*((x-c).^(a-1)).*((d-x).^(b-1));

 end

Would you kind test it first to see whether the code works and share with me every single line of the code? I am very grateful to your warm support. Notice also, c and d are restricted to 12-52, specially, c should be from 10-19, while d should between 35 and 52.

Kind regards,

Rdu end

[Rdu]

Hi Yu,

I also notice that you used the optimization of (err line) (i.e., least square) in myObj function file. I am not sure whether it can be used in fmincon. I have seen that most fmincon examples just provide function (foo) with extra code for minimization. The minimization codes (least square) are used in fminsearch. Am I right?

Furthermore, if I want to constrain for some parameters, say, the summation of a b c (just for an example, not for the present case) equal to 1. How can I add such a constraints in fmincon? Your expertise is highly appreciated.

Rdu

[Yu Jiang]

Hi Rdu

• Please notice that the function foo should have five arguments, namely, x,a,b,c,d, not x,p. However, the nonlinear function foo gives complex values even within the specified area. To get the code to work, I only take the real parts of the output of the function foo. You can find the code I fixed in the attached. I suggest you carefully check the feasible range of parameters such that only real numbers can be obtained from foo.
• Here are my answers to your second comment.

1) The reason that I minimize the least squares error is because this problem is a data fitting problem, in which one would want to find optimal parameters that can minimize the difference between the foo(x,a,b,c,d) and y. Simply minimizing foo would not provide the results you want.

2) If other linear equality constraints needs to be imposed, such as a+b+c+d = 1, the matrices Aeq and Beq can be set as Aeq = [1 1 1 1] and Beq = 1. Then, you may specify Aeq and Beq when using fmincon (See Documentation).

• On the other hand, I am pretty sure this is a standard data fitting problem that can be solved by using the Curve Fitting Toolbox (see documentation) , and I suggest you try it. In particular, this example explains how to perform data fitting with custom equations.

-Yu

[Rdu]

Hi Yu,

This is great! It works!!

I have tried to understand the underlying reason. I played with your code little bit. What I do not fully understand is that when I try to give x=[12:52] in all three functions explicitly (not just giving it in the main file [I modified foo((x(i)) in myObj and nlcon functions), it do not work. The reason why I am doing so is that I have seen in many cases that foo and myObj two functions can be combined (at least I am using in this way with fminsearch, for example "[err, yfit]=myObj(p,y)" ). With your rich expertise in the area, if possible, could you kindly tell me why my approach to combine these two functions in one function (one file) does not work?

Tons of thanks!

Rdu

[Yu Jiang]

I would be happy to help you resolve the issue. Could you post the code you revised? Also, please provide detailed steps for me to reproduce the issue.

[Rdu]

Dear Yu,

Thank you for your generous support!

I have figured out my problem about how to combine two functions together. See my codes in the attached. Now I have another question about the inconsistency by using fminsearch and fmincon with the same dataset and same initial values for parameter.

I am attaching you my m files plus a figure for comparing the outputs between two methods. It is my understanding that the fitted value from fminsearch and fmincon should be very close or identical since they use the same data with the same formula and the same initial values. Your feedback is highly appreciated on this discrepancy!

Rdu

[Rdu]

Hi Yu,

There is a typo in my main file. This one is the correct one. Thank you for your time.

Rdu

[Yu Jiang]

Hi Rdu

I believe fmincon and fminsearch use different algorithms to solve optimization problems. As you may have noticed, fmincon and fminsearch are used for bounded and unbounded optimization problems, respectively.

By testing your code on my machine, I noticed that fminsearch took 0.12 second and stopped because the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-08 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-06

However, fmincon took 0.81sec and stopped due to the objective function is non-decreasing in feasible directions, to within the selected value of the function tolerance, and constraints are satisfied to within the default value of the constraint tolerance.

I think this is a proof to show the two functions have quite different search algorithms. In fact, if you compare est1 ans est2, you may find they are not too different. I guess by setting appropriate optimset, you may bring them closer.

I understand that, in your figure, the differences between the curves can be easily observed. I think that is due to the high non-linearity of the function, which may greatly amplify the slight difference between est1 and est2.

-Yu

[Rdu]

Thank you for the explanation, Yu! I appreciate your time invested in my question in the past week. Good wishes to you! Let's close this question.

Rdu

Answer 2

[Rdu]

Thank you, Yu! I think my question is different from those in documentation and most at this website. My dataset has two variables x=[10:52}, and y=[0.006, 0.017, 0.026,..., 0.003, 0.002, 0.001]. I am going to estimate a mixture Beta distribution. function is f(x)={g(a+b)/[g(a)*g(b)]}*[(d-c)^(-a+b-1)]*[(x-c)^(a-1)]*[(d-x)^(b-1)]. a and b should be greater than 10 but less than 52. And the estimated f(x) should be greater than 0. How to program it using fmincon?

Appreciation.

Rdu

[Rajalekshmi kishhore]

can we have a separate function for linear constraint as was done for nonlinear constrint.if not how to mention sum linear constraint