function [x lambda] = eqp(G,c, A, b) %Equality Constrained Quadratic Problem Solver (eqp) : %Arash Mehraban, University of Colorado, Boulder % A' = QR where Q'Q = I (Q=n x n and R = n x m) % _____ _m_ % | | |\ | % | | |_\_|m % |_____| |___| % _____ _m_ % | | | |\r1| % A'= |y |z | |_\_|m % |__|__| |_0_|n-m % % A' = yr1 % Az = r1'y'z =0 %Inputs: % G = m x m matrix % c = m x 1 vector % A = m x n matrix % b = m x 1 vector % %Outputs: % x : solution of Quadratic: 1/2x'Gx + x'c subject to Ax = b % lambda: Lagrange Multiplier(s) %compute x = x0 + zu; [q r] = qr(A'); % find the QR factorization of A' [m n] = size(A); % find the dimensions of A r1 = r(1:m,:); % compute r1 from R in QR y = q(:, 1:m); % compute y from Q in QR z = q(:, m+1:n); % produce z matrix based on Q chol(z'*G*z); % error: reduced Hessian z'Gz is not positive definite x0 = y * inv(r1')* b; %compute x0 x = x0 - z * inv(z'*G*z) * z'*(G*x0 +c'); % compute x %compute lambda = inv(r1')* y'* (c + G*x) lambda = inv(r1)* y'* (c' + G*x); %compute lambda end