Re: How to Solve a Quadratic Programming problem?

*To*: mathgroup at smc.vnet.net*Subject*: [mg4216] Re: How to Solve a Quadratic Programming problem?*From*: rubin at msu.edu (Paul A. Rubin)*Date*: Tue, 18 Jun 1996 03:25:59 -0400*Organization*: Michigan State University*Sender*: owner-wri-mathgroup at wolfram.com

In article <4pisve$jk5 at dragonfly.wolfram.com>, daniele rizzi <drizzi at chiostro.univr.it> wrote: ->I'd like to solve a QP problem, in the following form: -> ->min translate(x)*H*x -> ->subject to translate(E)*x = E_f -> \sum_{i=1}^n x_i = 1 -> 0 <= x_i <= 1 \forall i \in (1, n) -> ->where : x is a n-term solution vector; -> E is a n-term constrain vector; -> E_f is a parameter of the problem; -> H is a n x n-term Symmetric Positive (semi)Definite Matrix. -> ->Mathematica is quite good at solving linear constrained (LP) problems -> (via Simplex-based routines) and general unconstrained instances, ->but what about quadratic objective function? Is there some "ready-made" ->routine or have I to write down my own code to tackle with that? -> There is no built-in routine for this in the current version of Mathematica. I don't know what is on tap for version 3.0. You might contact Jean-Christophe Culioli (culioli at cas.ensmp.fr); I believe he has working code to solve nonlinear programs with constraints. If you need to do it yourself, a quick-and-dirty approach is to relax the equation constraints into the objective as penalty terms. The problem is that as the errors get close to zero, higher powers of them get very small, so convergence is not great (i.e. the answer you get may not be very accurate in satisfying the constraints). An example: In[1]:= H := {{5., 3., -1.}, {3., 4., 0.}, {-1., 0., 10.}}; e = {1., 4., 2.}; ef = 3.1; u = {1., 1., 1.}; xvec = Array[ x, 3 ] (* vector of variables *) Out[1]= {x[1], x[2], x[3]} In[2]:= xlist = {#, 0.5, 0., 1.}& /@ xvec (* argument list for FindMinimum *) Out[2]= {{x[1], 0.5, 0., 1.}, {x[2], 0.5, 0., 1.}, {x[3], 0.5, 0., 1.}} In[3]:= f = (* objective function with quadratic penalties *) Simplify[ xvec . (H . xvec) + mu1 (e . xvec - ef)^2 + mu2 (u . xvec - 1)^2 ]; In[4]:= s = (* note the huge values for the penalty coefficients mu1, mu2 *) FindMinimum[ Evaluate[ f /. {mu1 -> 10^12, mu2 -> 10^12} ], Evaluate[ Sequence@@xlist ] ] Out[4]= {2.66847, {x[1] -> 0.114286, x[2] -> 0.607143, x[3] -> 0.278571}} In[5]:= xx = xvec /. s[[2]]; u . xx (* check conformance to constraints *) e . xx Out[5]= 1. 3.1 This can be a bit dicey in general, particularly if your optimum lies on the boundary of the feasible region (FindMinimum generates an error message if it thinks it is about to step outside the region). You can also use ConstrainedMin in a cutting plane scheme: In[6]:= const = { e . xvec == ef, u . xvec == 1 }; (* equation constraints *) cuts = {}; (* cutting planes *) vlist = Append[ xvec, z ]; (* variables - z = objective value *) In[7]:= s = ConstrainedMin[ z, Union[ const, cuts ], vlist ] Out[7]= {0, {x[1] -> 0.3, x[2] -> 0.7, x[3] -> 0, z -> 0}} In[8]:= x0 = xvec /. s[[2]]; (* compute the correct value of the quadratic *) z0 = x0 . (H . x0) Out[8]= 3.67 In[9]:= cuts = Append[ cuts, z >= Simplify[ z0 + 2 (H . x0) . (xvec - x0) ] ] Out[9]= {z >= -3.67 + 7.2 x[1] + 7.4 x[2] - 0.6 x[3]} In[10]:= s = ConstrainedMin[ z, Union[ const, cuts ], vlist ] (* try again *) Out[10]= {0.13, {x[1] -> 0, x[2] -> 0.55, x[3] -> 0.45, z -> 0.13}} In[11]:= x0 = xvec /. s[[2]]; z0 = x0 . (H . x0) Out[11]= 3.235 In[12]:= cuts = Append[ cuts, z >= Simplify[ z0 + 2 (H . x0) . (xvec - x0) ] ] Out[12]= {z >= -3.67 + 7.2 x[1] + 7.4 x[2] - 0.6 x[3], z >= -3.235 + 2.4 x[1] + 4.4 x[2] + 9. x[3]} And so forth. You can delete some of the older cuts (first in, first out should work) to keep the problem size compact, and limit the number of iterations. The x values will converge to the constrained optimum (assuming H is positive semidefinite), and the z value will converge to the optimal value of the quadratic objective. x will always be feasible (but suboptimal); z will be superoptimal (optimistic). Good luck. Paul ************************************************************************** * Paul A. Rubin Phone: (517) 432-3509 * * Department of Management Fax: (517) 432-1111 * * Eli Broad Graduate School of Management Net: RUBIN at MSU.EDU * * Michigan State University * * East Lansing, MI 48824-1122 (USA) * ************************************************************************** Mathematicians are like Frenchmen: whenever you say something to them, they translate it into their own language, and at once it is something entirely different. J. W. v. GOETHE ==== [MESSAGE SEPARATOR] ====