Thread Subject:

How to get max value of a function subject some constraints in

This is a math problem in my project.
If theta_N=[a b c d]', which is variable. theta_0=[a0 b0 c0 d0]',which is constant and known. PN is a 4X4 positive definite symetric matrix, which is constant and known.
If Y=T^2*(a-a0)+T*(b-b0)+(c-c0)+t*(d-d0).
Here, T and t are constant and known value. Y is scalar.
I need get the max value of |Y|, and corresponding theta_N

The constaint is
(theta_N-theta_0)'*(PN^-1)*(theta_N-theta_0)<=alpha
alpha is a constant and known value.

I asked this question in sci.math, someone told me I should use Lagrange multiplier method. But I don't known how to use this method in Matlab because I am not familar with Matlab.

BTW, the actual problem is even more complicated. In my project, T represents temperature and changes along with the time. t represents time, and the unit is second. At last I need get the max value of the sum of |Y1|+|Y2|+|Y3|+...+|Yn|, when the time is 24 hours, which makes n=24*3600. And I need get corresponding theta_N.

hhwolf76 <james.zhou76@gmail.com> wrote in message <10587277.969.1240498965778.JavaMail.jakarta@nitrogen.mathforum.org>...
> This is a math problem in my project.
> If theta_N=[a b c d]', which is variable. theta_0=[a0 b0 c0 d0]',which is constant and known. PN is a 4X4 positive definite symetric matrix, which is constant and known.
> If Y=T^2*(a-a0)+T*(b-b0)+(c-c0)+t*(d-d0).
> Here, T and t are constant and known value. Y is scalar.
> I need get the max value of |Y|, and corresponding theta_N
>
> The constaint is
> (theta_N-theta_0)'*(PN^-1)*(theta_N-theta_0)<=alpha
> alpha is a constant and known value.
>
> I asked this question in sci.math, someone told me I should use Lagrange multiplier method. But I don't known how to use this method in Matlab because I am not familar with Matlab.
>
> BTW, the actual problem is even more complicated. In my project, T represents temperature and changes along with the time. t represents time, and the unit is second. At last I need get the max value of the sum of |Y1|+|Y2|+|Y3|+...+|Yn|, when the time is 24 hours, which makes n=24*3600. And I need get corresponding theta_N.

  Define the column vectors x = theta_N-theta_0, and R = [T^2;T;1;t]. Then your problem is equivalent to finding x such that Y^2 = (x'*R)^2 = x'*(R*R')*x is maximized subject to x'*(PN^-1)*x <= alpha. This is an eigenvector problem.

 [U,S,V] = svd(PN*R*R');
 x = sqrt(alpha/(V(:,1)'*PN^-1*V(:,1)))*V(:,1);

  The reasoning is this. By the way x is defined it has the property that

 x'*(PN^-1)*x = alpha

Also

 PN*R*R'*V(:,1) = D(1,1)*V(:,1)

where D(1,1) is the largest eigenvalue and V(:,1) the corresponding eigenvector. Multiplying on the left by V(:,1)'*PN^-1 gives

 V(:,1)'*R*R'*V(:,1) = D(1,1)*V(:,1)'*PN^-1*V(:,1)

Thus

 Y^2 = x'*R*R'*x = D(1,1)*x'*PN^-1*x = D(1,1)*alpha

Therefore Y^2 takes its maximum value with this maximum eigenvalue, D(1,1).

  You can then find [a;b;c;d] as an offset of [a0,b0,c0,d0] by x.

  Your informant in sci.math was correct. This is indeed a problem in Lagrange multipliers and we have selected the largest multiplier by choosing the largest eigenvalue.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsqdk1$f0e$1@fred.mathworks.com>...
> ......
> PN*R*R'*V(:,1) = D(1,1)*V(:,1)
>
> where D(1,1) is the largest eigenvalue and V(:,1) the corresponding eigenvector. Multiplying on the left by V(:,1)'*PN^-1 gives
>
> V(:,1)'*R*R'*V(:,1) = D(1,1)*V(:,1)'*PN^-1*V(:,1)
>
> Thus
>
> Y^2 = x'*R*R'*x = D(1,1)*x'*PN^-1*x = D(1,1)*alpha
> ......

  I mistakenly used D for S in the above reasoning. Sorry. Replace each D(1,1) by S(1,1).

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsqdk1$f0e$1@fred.mathworks.com>...
> Define the column vectors x = theta_N-theta_0, and R = [T^2;T;1;t]. Then your problem is equivalent to finding x such that Y^2 = (x'*R)^2 = x'*(R*R')*x is maximized subject to x'*(PN^-1)*x <= alpha. This is an eigenvector problem.
>
> [U,S,V] = svd(PN*R*R');
> x = sqrt(alpha/(V(:,1)'*PN^-1*V(:,1)))*V(:,1);
>
> The reasoning is this. By the way x is defined it has the property that
>
> x'*(PN^-1)*x = alpha
>
> Also
>
> PN*R*R'*V(:,1) = D(1,1)*V(:,1)
>
> where D(1,1) is the largest eigenvalue and V(:,1) the corresponding eigenvector. Multiplying on the left by V(:,1)'*PN^-1 gives
>
> V(:,1)'*R*R'*V(:,1) = D(1,1)*V(:,1)'*PN^-1*V(:,1)
>
> Thus
>
> Y^2 = x'*R*R'*x = D(1,1)*x'*PN^-1*x = D(1,1)*alpha
>
> Therefore Y^2 takes its maximum value with this maximum eigenvalue, D(1,1).
>
> You can then find [a;b;c;d] as an offset of [a0,b0,c0,d0] by x.
>
> Your informant in sci.math was correct. This is indeed a problem in Lagrange multipliers and we have selected the largest multiplier by choosing the largest eigenvalue.
>
> Roger Stafford

  There is still an error in what I sent. I should have used the generalized 'eig' function rather than 'svd'. Below is the correction. Define x and R as before. Then do:

 PI = inv(PN); % PN is invertible
 [V,D] = eig(R*R',PI); % Generalized eigenvector function
 [t,k] = max(abs(diag(D))); % Find maximum eigenvalue
 x = real(sqrt(alpha/(V(:,k)'*PI*V(:,k)))*V(:,k));

  The reasoning now goes very much the same as before. Because of the way x is defined it must satisfy

 x'*PI*x = alpha

Also each eigenvector (including the k-th) from eig(R*R',PI) must satisfy

 R*R'*V(:,k) = D(k,k)*PI*V(:,k)

Multiplying on the left by V(:,k)' gives

 V(:,k)'*R*R'*V(:,k) = D(k,k)*V(:,k)'*PI*V(:,k)

Thus

 Y^2 = x'*R*R'*x = D(k,k)*x'*PI*x = D(k,k)*alpha

Therefore Y^2 takes its maximum value as this maximum eigenvalue, D(k,k)*alpha.

  As before, [a;b;c;d] is found as the offset of [a0,b0,c0,d0] by x.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsqdk1$f0e$1@fred.mathworks.com>...
> Define the column vectors x = theta_N-theta_0, and R = [T^2;T;1;t]. Then your problem is equivalent to finding x such that Y^2 = (x'*R)^2 = x'*(R*R')*x is maximized subject to x'*(PN^-1)*x <= alpha. This is an eigenvector problem.
> ..........
> You can then find [a;b;c;d] as an offset of [a0,b0,c0,d0] by x.
> ..........

  I hate to keeping throwing modifications at you but here is another way you can solve your problem that doesn't involve an explicit computation of an inverse. Since PN is positive definite and symmetric, it must possess an inverse, but computing that inverse directly is often considered undesirable. Instead do the following:

 [V,D] = eig(PN*R*R');
 [d,k] = max(abs(diag(D))); % Find maximum eigenvalue
 v = V(:,k);
 x = real(sqrt(d*alpha/(v'*R*R'*v))*v);

(The real part is taken in x here to remove any tiny imaginary component due to round-off error in 'eig'. Ideally this eigenvector for the largest eigenvalue should be real with all other eigenvalues exactly zero.)

  Reasoning: Since d and v are corresponding eigenvalue and eigenvector of PN*R*R', they must have the property

 PN*R*R'*v = d*v

Multiply both sides of this equation on the left by v'*inv(PN) (which we can do in principle) to get

 v'*R*R'*v = d*v'*inv(PN)*v

Then

 x'*inv(PN)*x = d*alpha/(v'*R*R'*v) * v'*inv(PN)*v
  = d*alpha/(v'*R*R'*v) * 1/d*v'*R*R'*v
  = alpha

and consequently x satisfies the required constraint.

  Then Y^2 = (R'*x)^2 = x'*R*R'*x will be the maximum possible value of Y^2 subject to the condition x'*inv(PN)*x <= alpha

Roger Stafford

Thanks for answering my question.
Here is the extension of my question. The constraint is still (theta_N-theta_0)'*(PN^-1)*(theta_N-theta_0)<=alpha
For specific T and t, I can solve the theta_N and corresponding Y by using your solution.
But if T and t keep changing, how can I get the max value of the sum of |Y1+Y2+...+Yn| and corresponding theta_N? n could be 24*3600. There are lots of data here.

First solution jumping to my head is to solve the every single theta_N and calculate corresponding |Y1+Y2+...+Yn|, comparing them and find the max one. But it need a lot of calculation. Is there some simple solution for this? Thank you very much!

hhwolf76 <james.zhou76@gmail.com> wrote in message <30651405.54061.1241467743222.JavaMail.jakarta@nitrogen.mathforum.org>...
> Thanks for answering my question.
> Here is the extension of my question. The constraint is still (theta_N-theta_0)'*(PN^-1)*(theta_N-theta_0)<=alpha
> For specific T and t, I can solve the theta_N and corresponding Y by using your solution.
> But if T and t keep changing, how can I get the max value of the sum of |Y1+Y2+...+Yn| and corresponding theta_N? n could be 24*3600. There are lots of data here.
>
> First solution jumping to my head is to solve the every single theta_N and calculate corresponding |Y1+Y2+...+Yn|, comparing them and find the max one. But it need a lot of calculation. Is there some simple solution for this? Thank you very much!
----------
  It is important to distinguish between maximizing |Y1+Y2+...+Yn| and maximizing |Y1|+|Y2|+|Y3|+...+|Yn|, the latter of which you mentioned in your first query in this thread. The solution for the first of these can well make some of the Y's negative in reaching its optimum solution, which would influence the selection of that optimum, whereas negative values for some of the Y's don't matter in the second summation.

  To maximize |Y1+Y2+...+Yn| all you need to do is add up all the R's (using my notation for R = [T^2;T;1;t]), to get a new R and then apply the most recent eigenvector method I mentioned to that R to solve the problem. (Presumably your PN remains the same for all the Y's.)

  I don't know how to maximize |Y1|+|Y2|+|Y3|+...+|Yn| which is an L1 metric problem. One could maximize Y1^2+Y2^2+Y3^2+...+Yn^2, which is an L2 problem involving a different method of calculating the final R.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gtnrq2$a1v$1@fred.mathworks.com>...
> ........
> I don't know how to maximize |Y1|+|Y2|+|Y3|+...+|Yn| which is an L1 metric problem. ........
----------
  When I said "I don't know how to maximize |Y1|+|Y2|+|Y3|+...+|Yn|", I should really have said that I don't know how to maximize it without using iterative techniques like those in the Optimization Toolbox with such functions as 'fmincon'. Your problem is a four-variable problem that 'fmincon' should be able to handle successfully by minimizing the negative of the above sum subject to your given constraints. Of course non-iterative techniques are usually preferable when available since they don't involve the uncertainties of making initial estimates for solutions and probably require fewer steps.

Roger Stafford

Sorry, It has been several months, but I suddenly jump out a question about that.

How do you get
Y^2 = x'*R*R'*x = D(1,1)*x'*PN^-1*x = D(1,1)*alpha

from

V(:,1)'*R*R'*V(:,1) = D(1,1)*V(:,1)'*PN^-1*V(:,1)

obviously, you have already imply
x = sqrt(alpha/(V(:,1)'*PN^-1*V(:,1)))*V(:,1);

But what if x is not?

If x = sqrt(alpha/(w'*PN^-1*w))*w;
w is random vector (not eigenvector), still
x'*(PN^-1)*x = alpha;

How could we know this
x'*R*R'*x <= D(1,1)*alpha

Thank you very much

Hui