Wikimization - New pages [en]

Linear Matrix Inequality II

Ranjelin — Sun, 15 Feb 2026 22:19:27 GMT

Summary: /* Convexity of the LMI constraint */

In convex optimization, a '''linear matrix inequality (LMI)''' is an expression of the form
: <math>LMI(y):=A_0+y_1A_1+y_2A_2+\ldots+y_m A_m\geq0\,</math>
where
* <math>y=[y_i\,,\,i\!=\!1\ldots m]</math> is a real vector,
* <math>A_0\,, A_1\,, A_2\,,\ldots\,A_m</math> are symmetric matrices in the subspace of <math>n\times n</math> symmetric matrices <math>S^n</math>,
* <math>B\geq0 </math> is a generalized inequality meaning <math>B</math> is a positive semidefinite matrix belonging to the positive semidefinite cone <math>S_+</math> in the subspace of symmetric matrices <math>S</math>.

This linear matrix inequality specifies a convex constraint on ''y''.

== Convexity of the LMI constraint ==
<math>LMI(y)\geq 0</math> is a convex constraint on ''y'' which means membership to a dual (convex) cone as we now explain: '''('''[http://meboo.convexoptimization.com/Meboo.html Dattorro, Example 2.13.5.1.1]''')'''

Consider a peculiar vertex-description for a [[Convex cones|convex cone]] defined over the positive semidefinite cone

'''('''instead of the more common nonnegative orthant, <math>x\geq0</math>''')''':

for <math>X\!\in S^n</math> given <math>\,A_j\!\in S^n</math>, <math>\,j\!=\!1\ldots m</math>

<math>\mathcal{K}=\left\{\left[\begin{array}{c}\langle A_1\,,\,X\rangle\\ . \\ . \\ . \\\langle A_m\;,\,X\rangle\end{array}\right]|\;X\succeq0\right\}\subseteq{\mathbb{R}}^m</math>

<math>=\left\{\left[\begin{array}{c}{\text svec}(A_1)^T\\ . \\ . \\ . \\{\text svec}(A_m)^T\end{array}\right]{\text svec}X\;|\;X\!\succeq_{\!}0\right\}</math>

<math>:=\;\{A\,{\text svec}X\;|\;X\!\succeq_{\!}0_{}\}</math>

where
*<math>A\!\in_{}\!R^{m\times n(n+1)/2}</math>,
*symmetric vectorization svec is a stacking of columns defined in '''('''[http://meboo.convexoptimization.com/Meboo.html Dattorro, ch.2.2.2.1]''')''',
*<math>A_0=\mathbf{0}</math> is assumed without loss of generality.

<math>\mathcal{K}</math> is a [[Convex cones|convex cone]] because

<math>A\,\textrm{svec}{X_{_{p_1}}}_{\,},_{_{}}A\,\textrm{svec}{X_{_{p_2}}}\!\in\mathcal{K}\;\Rightarrow\;
A(\zeta_{\,}\textrm{svec}{X_{_{p_1}\!}}+_{}\xi_{\,}\textrm{svec}{X_{_{p_2}}})\in_{}\mathcal{K}
\textrm{\;\;for\,all\;\,}\zeta_{\,},\xi\geq0</math>

since a nonnegatively weighted sum of positive semidefinite matrices must be positive semidefinite.

Now consider the (closed convex) dual cone:

<math>\begin{array}{rl}\mathcal{K}^*
\!\!\!&=_{}\left\{_{}y\;|\;\langle z\,,\,y_{}\rangle\geq_{}0\,\,\textrm{for\,all}\;\,z\!\in_{_{}\!}\mathcal{K}_{}\right\}\subseteq_{}R^m\\
&=_{}\left\{_{}y\;|\;\langle z\,,\,y_{}\rangle\geq_{}0\,\,\textrm{for\,all}\;\,z_{\!}=_{\!}A\,\textrm{svec}X\,,\,X\geq0_{}\right\}\\
&=_{}\left\{_{}y\;|\;\langle A\,\textrm{svec}X\,,\,y_{}\rangle\geq_{}0\,\,\textrm{for\,all}\;\,X\!\geq_{_{}\!}0_{}\right\}\\
&=\left\{y\;|\;\langle\textrm{svec}X\,,\,A^{T\!}y\rangle\geq_{}0\;\;\textrm{for\,all}\;\,X\!\geq_{\!}0\right\}\\
&=\left\{y\;|\;\textrm{svec}^{-1}(A^{T\!}y)\geq_{}0\right\}
\end{array}</math>

that follows from Fejer's dual generalized inequalities for the positive semidefinite cone:

* <math>Y\geq0\;\Leftrightarrow\;\langle Y\,,\,X\rangle\geq0\;\;\textrm{for\,all}\;\,X\geq0</math>

This leads directly to an equally peculiar halfspace-description

<math>\mathcal{K}^*=\{y\!\in R^m\;|\,\sum_{j=1}^my_jA_j\geq_{}0_{}\}</math>

The summation inequality with respect to the positive semidefinite cone
is known as a ''linear matrix inequality''.

== LMI Geometry ==

Although matrix <math>\,A\,</math> is finite-dimensional, <math>\mathcal{K}</math> is generally not a polyhedral cone
(unless <math>\,m\,</math> equals 1 or 2) simply because <math>\,X\!\in S_+^n\,.</math>

Relative interior of <math>\mathcal{K}</math> may always be expressed
<math>{\rm rel\,int}\,\mathcal{K}=\{A\,{\rm svec}X\;|\;X\!>_{\!}0_{}\}.</math>

Provided the <math>\,A_j</math> matrices are linearly independent, then
<math>{\rm rel\,int}\,\mathcal{K}={\rm int}\,\mathcal{K}</math>

meaning, cone <math>\mathcal{K}</math> interior is nonempty; implying, dual cone <math>\mathcal{K}^*</math> is pointed ([http://meboo.convexoptimization.com/Meboo.html Dattorro, ch.2]).

If matrix <math>\,A\,</math> has no nullspace, then
<math>\,A\,{\rm svec}X\,</math> is an isomorphism in <math>\,X\,</math> between the positive semidefinite cone <math>S_+^n</math> and range <math>\,\mathcal{R}(A)\,</math> of matrix <math>\,A.</math>

That is sufficient for [[Convex cones|convex cone]] <math>\,\mathcal{K}\,</math> to be closed, and necessary to have relative boundary
<math>{\rm rel}\,\partial^{}\mathcal{K}=\{A\,{\rm svec}X\;|\;X\!\geq_{\!}0\,,\,\det X=0_{}\}.</math>

<br>
Relative interior of the dual cone may always be expressed
<math>{\rm rel\,int}\,\mathcal{K}^*=\{y\!\in_{}\!R^m\;|\,\sum_{j=1}^my_jA_j>_{}0_{}\}.</math>

When the <math>A_j</math> matrices are linearly independent, function <math>\,g(y)_{\!}:=_{_{}\!}\sum y_jA_j\,</math> is a linear bijection on <math>R^m.</math>

Inverse image of the positive semidefinite cone under <math>\,g(y)\,</math>
must therefore have dimension equal to <math>\dim\!\left(\mathcal{R}(A^{\rm T})_{}\!\cap{\rm svec}\;S_+^{_{}n}\right)</math>

and relative boundary
<math>{\rm rel\,}\partial\mathcal{K}^*=\{y\in R^m \;|\,\sum_{j=1}^m y_j A_j\geq 0\,,\,\det\sum_{j=1}^m y_j A_j=0\}</math>

When this dimension is <math>\,m\,</math>, the dual cone interior is nonempty
<math>{\rm rel\,int}\,\mathcal{K}^*={\rm int}\,\mathcal{K}^*</math>

and closure of convex cone <math>\mathcal{K}</math> is pointed.

== Applications ==
There are efficient numerical methods to determine whether an LMI is feasible (''i.e.'', whether there exists a vector <math>y</math> such that <math>LMI(y)\geq0</math> ), or to solve a convex optimization problem with LMI constraints.
Many optimization problems in control theory, system identification, and signal processing can be formulated using LMIs. The prototypical primal and dual semidefinite program are optimizations of a real linear function respectively subject to the primal and dual [[Convex cones|convex cones]] governing this LMI.

== External links ==
* S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, [http://www.stanford.edu/~boyd/lmibook Linear Matrix Inequalities in System and Control Theory]

* C. Scherer and S. Weiland, [http://w3.ele.tue.nl/nl/cs/education/courses/hyconlmi Course on Linear Matrix Inequalities in Control], Dutch Institute of Systems and Control (DISC).

THD from Mapping Coefficients

Ranjelin — Wed, 22 Jan 2025 22:58:38 GMT

Summary: Matlab program for THD calculation

<pre>
%THD from mapping coefficients. Submeasurable Op Amp Distortion, section 3
function [thd dc] = thdxi(xi, E, precision);
mp.Digits(precision); %Advanpix MCT
Nh = numel(xi);
two = mp('2');
harmonic = zeros(Nh+1,1,'mp');
for n = 1:Nh
tn = xi(n)*E^n/two^n;
for ell = 0:n
tell = tn*nchoosek(mp(n), mp(ell));
idx = abs(n - 2*ell) + 1;
harmonic(idx) = harmonic(idx) + tell;
end
end
thd = sum(harmonic(3:end).^2)/harmonic(2)^2;
dc = harmonic(1);
end
</pre>

Diode Circuit Analysis

Ranjelin — Thu, 16 Jan 2025 23:51:32 GMT

Summary:

[[Image:Diode.jpg|850px]]
[[Image:Cyclic.jpg|850px]]
<pre>
function [vl Cdmax Cdmin Cdmed] = diodeMap(vi, Rs, Rl, verbose, fundfreq, Fs, precision) %from Convex Optimization & Euclidean Distance Geometry
mp.Digits(precision); %Advanpix MCT
k = 1.380649e-23; %boltzmann constant % https://www.nist.gov/si-redefinition/kelvin-boltzmann-constant
T = 298; %temperature kelvin %Neudeck p.5
q = 1.602176634e-19; %electron charge %NIST
is = 76.9e-12; %circuitlab diodeInduced
eta = 1.45; %circuitlab diodeInduced
tau = 4.32e-6; %transport time. max exponent (observed) e-4. Increasing capacitance observes more nonlinearity in vL FFT.
Rd = 0.042; %diode internal resistance
nVt = eta*k*T/q; %[V] thermal voltage for particular diode characterized by eta
isR = is*(Rs*Rl/(Rs + Rl) + Rd);

dvd = 1000*randn(size(vi)); %perturbation of derivative proves that iteration converges to same median capacitance
for i=1:5 %3 iterations will introduce randomness into mean impedance calculation
Lw = exp((isR + vi*Rl/(Rl + Rs))/nVt)*isR.*(1 + tau*dvd/nVt)/nVt;
vd = isR - nVt*Lambert_W(Lw) + vi*Rl/(Rl + Rs);
dvd = [vd(2)-vd(end); vd(3:end)-vd(1:end-2); vd(1)-vd(end-1)]/2; %approximate the derivative first time around, otherwise dvd becomes near 0. Assumes periodicity.
if verbose
Cd = tau*is*exp(vd/nVt)/nVt;
Cdmed = median(Cd);
disp([num2eng(Cdmed,false,false,false,16) ' farad']);
end
end
if verbose
Cdmax = max(Cd);
Cdmin = min(Cd);

figure(20); plot(0:numel(Cd)-1, Cd);
title('Cd over time'); xlabel('sample number'); ylabel('Farad'); set(get(gca,'ylabel'),'rotation',0);
drawnow
else
Cd = tau*is*exp(vd/nVt)/nVt;
end
iC = Cd.*dvd*Fs; %diffusion capacitor current. Normalization by Fs makes current independent of sample rate
id = is*(exp(vd/nVt) - 1); %ideal Shockley equation. This is already independent of Fs.
vl = (id + iC)*Rd + vd;

if verbose
figure(21); plot(0:numel(iC)-1, iC);
title('Cd current'); xlabel('sample number'); ylabel('Amp'); set(get(gca,'ylabel'),'rotation',0);
figure(22); plot(0:numel(id)-1, id);
title('diode current'); xlabel('sample number'); ylabel('Amp'); set(get(gca,'ylabel'),'rotation',0);
idx = find(iC ~= 0);
disp(['mean capacitor impedance = ' num2str(orosumvec(vd(idx)./iC(idx),1)/numel(idx),'%3.4e')])
% disp(['mean capacitor |impedance| = ' num2str(orosumvec(abs(vd(idx)./iC(idx)),1)/numel(idx),'%3.4e')])
t = median(vd./iC);
if t < 0, spc=[]; else spc=' '; end
disp(['median capacitor impedance =' spc num2str(t)]) %median impedance is always far smaller than mean despite mp.Digits precision
medcapimp = median(abs(vd./iC));
disp(['median capacitor |impedance| = ' num2str(medcapimp,'%3.4e')])
disp(['empirical constant c = ' num2str(medcapimp*(2*pi*fundfreq*Cdmed),'%3.15e')])

figure(23);
histogram(Cd,50,'Normalization','probability'); %for book illustration
% histogram(Cd,round(3.25*48000)); %for counting smallest diffusion capacitors
title('histogram C_d'); xlabel('bin'); ylabel('relative count'); %set(get(gca,'ylabel'),'rotation',0);
drawnow
end
end

% https://www.mathworks.com/matlabcentral/fileexchange/43419-the-lambert-w-function
function w = Lambert_W(x,branch)
% Lambert_W Functional inverse of x = w*exp(w).
% w = Lambert_W(x), same as Lambert_W(x,0)
% w = Lambert_W(x,0) Primary or upper branch, W_0(x)
% w = Lambert_W(x,-1) Lower branch, W_{-1}(x)
%
% See: http://blogs.mathworks.com/cleve/2013/09/02/the-lambert-w-function/

% Effective starting guess
if nargin < 2 || branch ~= -1
w = ones(size(x)); % Start above -1
else
w = -2*ones(size(x)); % Start below -1
end
v = inf*w;

% Halley's method
c = 0;
limit = 100;
while any(abs(w - v)./abs(w) > 1e-15) && c < limit %changed magic tolerance. Was 1e-8. -JonD
v = w;
e = exp(w);
f = w.*e - x; % Iterate, to make this quantity zero
w = w - f./((e.*(w+1) - (w+2).*f./(2*w + 2)));
c = c + 1;
end
if c >= limit
disp('%%%%%%%%%%%%%%%%%%%%%% Warning: Lambert limit reached %%%%%%%%%%%%%%%%%%%%%')
end
end
</pre>