How to Solve Linear Complementarity Problem (LCP) with MATLAB ?

I want to use MATLAB for solving Linear Complementarity Problem (LCP) as defined in the following:

what optimization function (e.g., fmincon,lsqnonneg,etc.) should I use, or someone could give me some example.

[Matt J]

You don't appear to have an objective function. Are you just trying to find an arbitrary (ω, z) satisfying these (in)equalities?

Answer 1

[Matt J]

You could do it with quadprog, but you may as well just deploy it as an optimproblem.

n=5; q=-rand(n,1); M=rand(n); M=(M+M')/2;

smallTolerance=1e-6;

z=optimvar('z',n,'Lower',0);

prob=optimproblem('Objective', z'*M*z + z'*q, "Constraints", M*z+q>=0);

[sol,fval]=solve(prob);

Solving problem using quadprog.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

fval

fval = 7.2164e-16

assert(fval<=smallTolerance, "No solution found")

eig(M)', %Not PSD

ans = 1×5

-0.3846 -0.2232 0.0577 0.4325 2.3596

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

zOptimal=sol.z    %Solutions

zOptimal = 5×1

0.0000 0.8317 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

omegaOptimal=M*zOptimal+q

omegaOptimal = 5×1

0.1923 0 0.1923 0.3776 0.6387

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Torsten]

Which solver will be invoked by "solve" ?

Will "quadprog" work if M is not positive-semidefinite ?

Trying

https://uk.mathworks.com/matlabcentral/fileexchange/41485-a-pivoting-algorithm-solving-linear-complementarity-problems

might be an alternative.

[Matt J]

Quadprog will "work", even if M is not psd, but of couse it might then not find the global solution. I've modified my answer with an example where it succeeds, in spite of this.

[Torsten]

Don't you have to use

n=5; q=-rand(n,1); M=rand(n); Msym=(M+M')/2;
smallTolerance=1e-6;
z=optimvar('z',n,'Lower',0);
prob=optimproblem('Objective', z'*Msym*z + z'*q, "Constraints", M*z+q>=0);
[sol,fval]=solve(prob);
fval
assert(fval<=smallTolerance, "No solution found")
eig(Msym)', %Not PSD
zOptimal=sol.z    %Solutions
omegaOptimal=M*zOptimal+q

for M being a general nxn matrix ?

[Chuguang Pan]

@Matt J. Thanks for your detailed reply, your answer is helpful.