Solve an optimal control problem with the Wolfram Language?

I have Mathematica version 8, but will gladly get the latest version if it is necessary to accomplish the following goal: find two control functions that maximize a certain functional in the control functions and 2 other functions subject to inequality constraints on the control functions and also non-linear differential equation constraints involving all 4 functions. If this is not possible, I will settle for solving for only one control function. Can someone please give me an idea of what must be done?

I have been able to do some constrained functional optimizations in Mathematica by using time-dependent Lagrange multipliers, as described in Courant & Hilbert, Methods of Mathematical Physics, in Chapter, Section 7, Variational Problems With Subsidiary Conditions. Unfortunately, I cannot at this time locate one of those calculations to post here.

I don’t currently have Courant and Hilbert but I’m guessing that they are talking of a way to ensure the differential equation constraints. Of course the Lagrange multipliers are time-dependent. They are variables which have their own Hamilton Jacobi equation, which is just the dif eq constraints, and then they appear as dependent variables in the other H-J equations. The thing that concerns me is not the differential equation constraints (which can be handled with Lagrange multipliers), but the inequality constraints. I suppose I can introduce some auxiliary function into the Lagrangian which gives a high penalty for going outside the domain allowed by the inequality?

Inequality constraints can be converted to equality constraints with the addition of slack or surplus terms, which are generally the square of a new variable, to insure they stay positive.

I don’t understand how a slack term leads to an inequality. For example if I want variable y<1 then what term would you add to Lagrangian?

expr >= 0 should be changed to expr - s^2 == 0. expr <= 0 should be changed to expr + s^2 == 0

And s is a dynamic variable? So a new term in the Lagrangian would be lambda (y-1+s^2) for satisfying y<1?

I think I see: s is just a change of variable from y, not a dynamic variable with a separate Hamilton Jacobi equation.

This doesn’t work if I have another inequality to satisfy, such as y>0. In that case we can change the variable to y=(1+Tanh[s])/2 or s=InvTanh[2y-1)

There are some demonstration on this topic on the Wolfram Demonstrations Project, e.g.
Moon Landing Simulation. I would think that they may give you an insight on your problem.