I need an inverse in my model

The inverse operator is not supported in the YALMIP language as it is an extremely complex operation which does not have a nice representation. Unless you can reformulate your model to completely eliminate the need for the inverse (typically Schur complements, congruence transformations or variable changes) you will have to accept that you have to convert your nasty inverse to a (slightly less) nasty alternative constraint.

Bilinear reformulation

The standard alternative to a model involving \(X^{-1}y\) is to replace this expression with a new variable \(z\) and add the bilinear constraint \(Xz = y\) to your model.

Note that invertability is lost in this form though (for instance \((X=0, y=0, z=0\) is feasible). Unfortunately, adding a non-singularity constraint is even worse than trying to use an inverse.

One approach to keep non-singularity is to explicitly introduce a matrix \(Y\) representing the inverse, and add the constraint \(XY = I\) together with \(Xz == y\) or \(z == Yy\).

On the topic of non-singular matrices, it is worth mentioning that the set of singular matrices is of measure 0, which means that any singular matrix can be perturbed with an infinitely small perturbation to render it non-singular, and since we work with solvers with finite numerical precision (typically in the order of \(10^{-8}\) or so) it would make no sense to try to add a constraint that a matrix is non-singular. Your model somehow has to implictly encode solutions which are significantly non-singular by construction, if this is important to you.

MILP reformulation

A reformulation which can work if the matrix only depends on a small number of binary variables is to perform an explicit enumeration of all possible combinations. The model for \(z = X^{-1}y\) would then be

% Find all involved variables (assumed binary)
involved = recover(depends(X));
n = length(involved)

% Loop through all combinatorial cases and define relationship by implication
cases = binvar(1,2^n-1);
Model = [sum(cases) == 1]
for i = 0:2^n-1
    xi = dec2decbin(i,n)';
    assign(involved,xi);
    V = value(X);
    if rank(V) == n
        Model = [Model, implies(cases(i), [z == inv(V)*y, x == xi])];
    else
        Model = [Model, cases(i) == 0];
    end
end

Blackbox approach

If you only need some scalar function which involves an inverse of a matrix which depends on some variable, you can use a blackbox approach. The drawback is of course that you are hiding structure completely, and your arrive at a general nonlinear nonconvex nonsmooth program.

q = [2;1];
f = @(x)(q'*inv([1 + sin(1+x(1))^2 2;2 4 + cos(x(1)-x(2))^2])*q);

x = sdpvar(2,1);
y = blackbox(f,x);

optimize([],y);
value(y)
f(value(x))