How to use matlab in order to solve an optimization problem with integral and limit constrains?

Hello,

So, given the equation below, does anyone know how to introduce the constrains in Matlab in order to solve for the constants B and P?

Y(x) = (T(X) - B)*(X/2*pi)^4-P; %Equation

lim Y(X) = 0; % First constrain x --> infinity

Integral Y(X) dX = 0 % Second constrain X=0 to infinity

[Walter Roberson]

What is known about T?

[Manuel Oliveira]

T is an arbitrary function of X

[John D'Errico]

Is T KNOWN, but arbitrary? Anyway, if T is unknown, just some general function, then it will be impossible to solve for B and P, as functionals of T.

It also seems that the behavior of T is limited by your constraints. For example, what is the limit of T as x--> inf?

Answer 1

[Walter Roberson]

There are multiple solutions I think. One of them looks to be P=0 and B=0 and T(x) is 1/X^N where N is odd and at least 5

[Walter Roberson]

"I want to know how introduce in the optimization toolbox the equation with the constrains (at least the integral) so that it can minimize and calculate the values of B and P."

guess = rand(1,2) * 100;
Y = @(X, B, P) (T(X) - B) .* (X/2*pi).^4 - P;
IY = @(B, P) integral( @(X) Y(X, B, P), 0, inf);
fmincon( @(BP) IY(BP(1), BP(2)).^2, guess)

whether lim(Y) goes to 0 at X going to infinity depends upon the form of T rather than the values of B and P.

Minimizing the square of the integral is equivalent to trying to get the integral to be 0.

Answer 2

[John D'Errico]

This really is not a MATLAB question. Simple analysis and reasoning seems to suffice.

Assuming you have some KNOWN function T(x), then consider what happens at infinity for Y. We are told that Y(X)must approach zero.

Y(x) = (T(X) - B)*(X/2*pi)^4-P

In order for this to be zero as X-->inf, we need to know what is the limit of T(X) as x--> inf. Call it T_inf.

If T_inf is +/- inf, then nothing you put in there for B or P will cancel out that product of infs.

So now look at the product

(T(X) - B)*(X/2*pi)^4

If this is to be finite or zero, then we need to think about how T behaves as X-->inf.

Suppose that T(X) behaves like k*X^(-p)+T_inf. The important case is p=4. If T(X)-B goes to zero more slowly than X^4 goes to inf, then the product will diverge. No constants for B or P can exist to let Y approach zero. If T(X)-B goes to zero more rapidly than X^4, then P will be zero. In that case, it is also clear that B must be simply the limit of T as X-->inf. So we would have B=T_inf.

If T acts like k*X^(-4)+T_inf at inf, then again, we must have B=T_inf, but only in this case will P will be non-zero.

p = -k/2*pi

So there is no need to estimate B OR P. Both are given by the behavior of T(X) as X approaches inf.

In fact, the second integral constraint is not even needed. It simply serves to over-constrain the parameters.

If the behavior of T is completely unknown, then we cannot infer values for B or P, since they are completely determined by the behavior of T as x-->inf.

[Walter Roberson]

In order for the integral to be 0, it has to have regions of opposite sign, and the behavior at -inf has to mirror the behaviour at +inf -- since Y(x) goes to 0 at +inf, it must also go to 0 at -inf but with opposite sign. The (X/2*pi)^4 part goes to infinity in both directions so (T(X) - B) must go to 0 in both directions and must go faster that (X/2*pi)^4 goes to 0. With the k*X^(-p)+T_inf model, that requires p be odd and at least 5.

If p = 4 then X^(-p) is always positive. Hmmm, maybe you could supply the negative region by virtue of the -B, so you have a positive region near 0. Ah, but integral of k*x^(-4) crossing 0 is infinity * sign(k) so No, that cannot work.

[Walter Roberson]

I was mis-reading, I thought it was integral from -inf to +inf

Answer 3

[Manuel Oliveira]

So, this is what I came up with. However I need to stop the minimization of the integral when the value reaches 0 (or a very close value to 0).

Any help?

x0 = [0.01,0.01];
fun = @(y)trapz(x,((y(2)-(T-y(1)).*(x./(2*pi)).^4)));
      B=y(1);
      P=y(2);
problem = createOptimProblem('fmincon','objective',fun,...
          'x0',x0);
gs = GlobalSearch;
[x,fg,flg,og] = run(gs,problem);

[Walter Roberson]

What is your x in fun?

Why are you using "x0" as a name for the initial guesses for y ?

In your original question, T was a constant, but here you have used it as either a scalar or a vector the same length as x. Is what you write here as T equivalent to your original T evaluated at x, which I would find clearer expressed as Tx (especially if you had just assigned that in a call just before you created fun)

[Walter Roberson]

This version of the problem does not express your two required constraints, and so is minimized when P approaches -infinity provided that no element of T is infinity (because that would lead to NaN.)