Sequential Quadratic Programming: Primal

Min (1/2)x.G.x+c.x satisfying (Subscript[a, i]).x=(Subscript[b, i]) for i∈ℰ  and (Subscript[a, i]).x≥(Subscript[b, i]) for i∈ℐ. 
where G is symmetric

Code

Test Matrices etc.

Clear[RandomSym, RandomVec, RandomMat]
RandomSym[n_, Distribution_:UniformDistribution[{1,2}]]:=Module[
{Q,dia},
Dia=SparseArray[Band[{1,1}]:>RandomReal[Distribution],{n,n}];
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]];
Q.Dia.Transpose[Q]
]
RandomVec[n_,Distribution_:UniformDistribution[{-1,1}]]:= RandomReal[Distribution,n]
RandomMat[{m_,n_},Distribution_:UniformDistribution[{-1,1}]]:= RandomReal[Distribution,{m,n}]

n=1200;
{
Histogram[Eigenvalues[RandomSym[n,NormalDistribution[2,0.1]]]],
Histogram[Eigenvalues[RandomSym[n]]]
}

KKTQuadratic Solver

The following takes in matrices and vectors for the Quadratic program and solves the linear system of KKT equations.  It assumes that A has full rank and that G is symmetric.  It gives the only KKT point if G is positive definate

KKTSolve[{G_,c_},{A_,b_}]:= Module[
{KKTMat,KKTRHS,xλStar,n=Length[c]},
KKTMat=ArrayFlatten[{{G, -Transpose[A]},{A,0}}];
KKTRHS=Join[-c, b];
xλStar=LinearSolve[KKTMat, KKTRHS];
{Take[xλStar,n],Drop[xλStar,n]}
]

First a 2D test.

G= {
{1.0,0},
{0,3.0}
};
c= {1.0,2.0};
A= {{1.0,1.0}};
b = {1.2};
{xStar,λStar}=KKTSolve[{G,c},{A,b}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -2,2},
{x2,-2,2},
ColorFunction->Hue,
Contours->120
],
ContourPlot[A.{x1,x2}-b==0,
{x1, -2,2},
{x2,-2,2},
ContourStyle->Thickness[0.01]
],
Graphics[{Red, PointSize[0.02],Point[xStar]}]
]

Now a 3D test with 1 constraint.

G= {
{1.0,0.01,0.3},
{0.01,3.0,0.1},
{0.1, 0.1,4.1}
};
c= {1.0,2.0,3.0};
A= {{1.0,1.0,0.2}};
b = {1.2};
{xStar,λStar}=KKTSolve[{G,c},{A,b}]
Show[
ContourPlot3D[
(1/2){x1,x2,x3}.G.{x1,x2,x3} + c.{x1,x2,x3},
{x1, -2,2},
{x2,-2,2},
{x3,-2,2},
ColorFunction->Hue,
Contours->23,
RegionFunction->Function[{x1,x2,x3,z},(A.{x1,x2,x3}-b)[[1]] >0]
],
Graphics3D[{Red, PointSize[0.02],Point[xStar]}],
AxesLabel->{"x1","X2","x3"}
]

Now a 3D test with 2 constraints.

G= {
{1.0,0.01,0.3},
{0.01,3.0,0.1},
{0.1, 0.1,4.1}
};
c= {1.0,2.0,3.0};
A= {{1.0,1.0,0.2},{-2.0,1.2,0.1}};
b = {1.2,0.2};
{xStar,λStar}=KKTSolve[{G,c},{A,b}]
Show[
ContourPlot3D[
(1/2){x1,x2,x3}.G.{x1,x2,x3} + c.{x1,x2,x3},
{x1, -2,2},
{x2,-2,2},
{x3,-2,2},
ColorFunction->Hue,
Contours->23,
RegionFunction->Function[{x1,x2,x3,z},(A.{x1,x2,x3}-b)[[1]] >0]
],
Graphics3D[{Red, PointSize[0.02],Point[xStar]}],
AxesLabel->{"x1","X2","x3"},
PlotLabel->"1st Constraint"
]
Show[
ContourPlot3D[
(1/2){x1,x2,x3}.G.{x1,x2,x3} + c.{x1,x2,x3},
{x1, -2,2},
{x2,-2,2},
{x3,-2,2},
ColorFunction->Hue,
Contours->23,
RegionFunction->Function[{x1,x2,x3,z},(A.{x1,x2,x3}-b)[[2]] >0]
],
Graphics3D[{Red, PointSize[0.02],Point[xStar]}],
AxesLabel->{"x1","X2","x3"},
PlotLabel->"2nd Constraint"
]
Show[
ContourPlot3D[
(1/2){x1,x2,x3}.G.{x1,x2,x3} + c.{x1,x2,x3},
{x1, -2,2},
{x2,-2,2},
{x3,-2,2},
ColorFunction->Hue,
Contours->23,
RegionFunction->Function[{x1,x2,x3,z},And[(A.{x1,x2,x3}-b)[[1]] >0,(A.{x1,x2,x3}-b)[[2]] >0]]
],
Graphics3D[{Red, PointSize[0.02],Point[xStar]}],
AxesLabel->{"x1","X2","x3"},
PlotLabel->"Both Constraints"
]

Looks to me like it is working for my examples.

Just as a check lets see what happens with a non PD G.  Of course it does not give the min but …

G= {
{-1.0,0},
{0,3.0}
};
c= {1.0,2.0};
A= {{1.0,1.0}};
b = {1.2};
{xStar,λStar}=KKTSolve[{G,c},{A,b}]
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2} + c.{x1,x2},
{x1, -3,3},
{x2,-3,3},
ColorFunction->Hue,
Contours->120
],
ContourPlot[A.{x1,x2}-b==0,
{x1, -3,3},
{x2,-3,3},
ContourStyle->Thickness[0.01]
],
Graphics[{Red, PointSize[0.02],Point[xStar]}]
]

Inequality Constrained problems: SQP

Here is a picture of a QP in 2D with 3 inequality constraints.

n=2;
G =RandomSym[n];
c=RandomVec[n];
mI=3;
{AI,bI}= {RandomMat[{mI,n}],RandomVec[mI]};
M=5;
Show[
ContourPlot[
(1/2){x1,x2}.G.{x1,x2}+ c.{x1,x2},
{x1,-M,M},{x2,-M,M},
Contours->60 ,
ColorFunction->ColorData["TemperatureMap"],
RegionFunction->Function[{x1,x2,f},Apply[And,Table[AI[[i]].{x1,x2}≥bI[[i]],{i,1,mI}]]]
],
Table[
ContourPlot[
AI[[i]].{x1,x2}-bI[[i]]==0,
{x1,-M,M},{x2,-M,M},
ContourStyle->{Thickness[0.01],Hue[i/mI]},
ContourShading->False
],
{i,mI}]
]

Here is a bad attempt in 3D showing the internal structure.

n=3;
G =RandomSym[n];
c=RandomVec[n];
mI=4;
{AI,bI}= {RandomMat[{mI,n}],RandomVec[mI]};
M=5;
Show[
ContourPlot3D[
(1/2){x1,x2,x3}.G.{x1,x2,x3}+ c.{x1,x2,x3},
{x1,-M,M},{x2,-M,M},{x3,-M,M},
Contours->24 ,
ColorFunction->ColorData["TemperatureMap"],
RegionFunction->Function[{x1,x2,x3,f},Apply[And,Table[AI[[i]].{x1,x2,x3}≥bI[[i]],{i,1,mI}]]]
],
ContourPlot3D[
Evaluate[Table[AI[[i]].{x1,x2,x3}-bI[[i]]==0,{i,mI}]],
{x1,-M,M},{x2,-M,M},{x3,-M,M},
Mesh->True,
MeshStyle->Opacity[0.2],
ContourStyle->Opacity[0.5]
]
]

Here is the same plot shoing only the external structure.

Show[
RegionPlot3D[
Apply[And,Table[AI[[i]].{x1,x2,x3}-bI[[i]]>=0,{i,1,mI}]],
{x1,-M,M},{x2,-M,M},{x3,-M,M},
Mesh->False,
PerformanceGoal->"Quality",
PlotPoints->80,
ColorFunction->Function[{x1,x2,x3},ColorData["TemperatureMap"][(1/2){x1,x2,x3}.G.{x1,x2,x3}+ c.{x1,x2,x3}]]
],
ContourPlot3D[
Evaluate[Table[AI[[i]].{x1,x2,x3}-bI[[i]]==0,{i,mI}]],
{x1,-M,M},{x2,-M,M},{x3,-M,M},
Mesh->True,
MeshStyle->Opacity[0.2],
ContourStyle->Opacity[0.5]
]
]