Extending the Law of Cosines formula to quadrilaterals (and even polygon?)
In quadrilateral $ABCD$, in addition to $a = BC$, $b = CD$, $c = AB$, $d = AD$, also set: $$ AG=e,\quad CG=f,\quad BG=g,\quad DG=h,\quad \angle BGC=\alpha, $$ where $G$ is the intersection point of diagonals $AC$ and $BD$. By the cosine law we get then: $$ \begin{align} a^2 &=f^2+g^2-2fg\cos\alpha \\ d^2 &=e^2+h^2-2eh\cos\alpha \\ b^2 &=f^2+h^2+2fh\cos\alpha \\ c^2 &=e^2+g^2+2eg\cos\alpha \\ \end{align} $$ and from that we obtain: $$ b^2+c^2-a^2-d^2 = 2(fh+eg+fg+eh)\cos\alpha=2(e+f)(g+h)\cos\alpha, $$ which is precisely your formula.
Let us consider convex $n$-gon $A_1A_2\cdots A_n$ where $\overline{A_jA_{j+1}}=a_j$ with $\angle{A_jA_{j+1}A_{j+2}}=\theta_j$.
Now, let us put our $n$-gon on the $xy$ plane in the following way :
$A_1$ is at the origin
The side $A_1A_2$ is on the $x$-axis
The $x$-coordinate of $A_2$ is positive
The $y$-coordinate of $A_3$ is positive.
Here, if we consider the projection of each side on the $x$-axis, then we get $$a_1+a_2\cos(\pi-\theta_1)+a_3\cos(2\pi-(\theta_1+\theta_2))+\cdots +a_n\cos((n-1)\pi-(\theta_1+\theta_2+\cdots +\theta_{n-1}))=0$$ which can be written as $$a_1=\sum_{k=1}^{n-1}(-1)^{k+1}a_{k+1}\cos\bigg(\sum_{j=1}^{k}\theta_j\bigg)\tag1$$
Similarly, if we consider the projection of each side on the $y$-axis, then we get $$a_2\sin(\pi-\theta_1)+a_3\sin(2\pi-(\theta_1+\theta_2))+\cdots +a_n\sin((n-1)\pi-(\theta_1+\theta_2+\cdots +\theta_{n-1}))=0$$ which can be written as $$0=\sum_{k=1}^{n-1}(-1)^{k+1}a_{k+1}\sin\bigg(\sum_{j=1}^{k}\theta_j\bigg)\tag2$$
From $(1)(2)$, we obtain $$a_1^2+0^2=\bigg(\sum_{k=1}^{n-1}(-1)^{k+1}a_{k+1}\cos\bigg(\sum_{j=1}^{k}\theta_j\bigg)\bigg)^2+\bigg(\sum_{k=1}^{n-1}(-1)^{k+1}a_{k+1}\sin\bigg(\sum_{j=1}^{k}\theta_j\bigg)\bigg)^2$$ which can be written as $$a_1^2=\sum_{k=1}^{n-1}a_{k+1}^2+\sum_{1\le p\lt q\le n-1}\bigg(2(-1)^{p+1}a_{p+1}\cos\bigg(\sum_{j=1}^{p}\theta_j\bigg)\times (-1)^{q+1}a_{q+1}\cos\bigg(\sum_{j=1}^{q}\theta_j\bigg)+2(-1)^{p+1}a_{p+1}\sin\bigg(\sum_{j=1}^{p}\theta_j\bigg)\times (-1)^{q+1}a_{q+1}\sin\bigg(\sum_{j=1}^{q}\theta_j\bigg)\bigg)$$ i.e. $$a_1^2=\sum_{k=1}^{n-1}a_{k+1}^2+\sum_{1\le p\lt q\le n-1}2(-1)^{p+q}a_{p+1}a_{q+1}\bigg(\cos\bigg(\sum_{j=1}^{p}\theta_j\bigg)\cos\bigg(\sum_{j=1}^{q}\theta_j\bigg)+\sin\bigg(\sum_{j=1}^{p}\theta_j\bigg)\sin\bigg(\sum_{j=1}^{q}\theta_j\bigg)\bigg)$$ i.e. $$a_1^2=\sum_{k=1}^{n-1}a_{k+1}^2+\sum_{1\le p\lt q\le n-1}2(-1)^{p+q}a_{p+1}a_{q+1}\cos\bigg(\sum_{j=1}^{q}\theta_j-\sum_{j=1}^{p}\theta_j\bigg)$$ Therefore, we get $$\color{red}{a_1^2=\sum_{k=1}^{n-1}a_{k+1}^2+\sum_{1\le p\lt q\le n-1}2(-1)^{p+q}a_{p+1}a_{q+1}\cos\bigg(\sum_{j=p+1}^{q}\theta_j\bigg)}$$
For example, for pentagon $A_1A_2A_3A_4A_5\ (n=5)$, we get
$$\color{red}{a_1^2=a_2^2+a_3^2+a_4^2+a_5^2-2a_{2}a_{3}\cos(\theta_2)+2a_{2}a_{4}\cos(\theta_2+\theta_3)-2a_{2}a_{5}\cos(\theta_2+\theta_3+\theta_4)-2a_{3}a_{4}\cos(\theta_3)+2a_{3}a_{5}\cos(\theta_3+\theta_4)-2a_{4}a_{5}\cos(\theta_4)}$$
Added : One can get several formulas.
For quadrilateral $A_1A_2A_3A_4\ (n=4)$ :
If we change $(1)(2)$ to $$(1)\implies a_4\cos(\theta_4)=a_1-a_2\cos(\theta_1)+a_3\cos(\theta_1+\theta_2)$$$$(2)\implies a_4\sin(\theta_4)=a_2\sin(\theta_1)-a_3\sin(\theta_1+\theta_2)$$squaring and adding give $$a_4^2=a_1^2+a_2^2+a_3^2-2a_1a_2\cos(\theta_1)-2a_2a_3\cos(\theta_2)+2a_1a_3\cos(\theta_1+\theta_2)$$
If we change $(1)(2)$ to $$(1)\implies a_1+a_3\cos(\theta_1+\theta_2)=a_2\cos(\theta_1)+a_4\cos(\theta_4)$$$$(2)\implies a_3\sin(\theta_1+\theta_2)=a_2\sin(\theta_1)-a_4\sin(\theta_4)$$squaring and adding give $$a_1^2+a_3^2+2a_1a_3\cos(\theta_1+\theta_2)=a_2^2+a_4^2+2a_2a_4\cos(\theta_1+\theta_4)$$
If we change $(1)(2)$ to $$(1)\implies a_3\cos(\theta_1+\theta_2)-a_2\cos(\theta_1)=a_4\cos(\theta_4)-a_1$$$$(2)\implies a_3\sin(\theta_1+\theta_2)-a_2\sin(\theta_1)=-a_4\sin(\theta_4)$$squaring and adding give $$a_2^2+a_3^2-2a_2a_3\cos(\theta_2)=a_1^2+a_4^2-2a_1a_4\cos(\theta_4)$$
For pentagon $A_1A_2A_3A_4A_5\ (n=5)$ :
- If we change $(1)(2)$ to $$(1)\implies a_1-a_2\cos(\theta_1)+a_3\cos(\theta_1+\theta_2)=a_5\cos(\theta_5)-a_4\cos(\theta_4+\theta_5)$$$$(2)\implies a_2\sin(\theta_1)-a_3\sin(\theta_1+\theta_2)=a_5\sin(\theta_5)-a_4\sin(\theta_4+\theta_5)$$squaring and adding give $$a_1^2+a_2^2+a_3^2-2a_1a_2\cos(\theta_1)-2a_2a_3\cos(\theta_2)+2a_1a_3\cos(\theta_1+\theta_2)=a_4^2+a_5^2-2a_4a_5\cos(\theta_4)$$
Let $\vec{BC}=\vec{a},$ $\vec{CD}=\vec{b},$ $\vec{DA}=\vec{d}$ and $\vec{AB}=\vec{c}.$
Thus, since $$\vec{a}+\vec{c}=-\vec{b}-\vec{d},$$ we obtain: $$(\vec{a}+\vec{c})^2=(\vec{b}+\vec{d})^2,$$ which gives $$\vec{a}\vec{c}-\vec{b}\vec{d}=\frac{1}{2}(b^2+d^2-a^2-c^2).$$ In another hand, $$BD\cdot AC\cos\measuredangle BGC=\vec{DB}\cdot\vec{AC}=(\vec{c}+\vec{d})(\vec{c}+\vec{a})=$$ $$=c^2+\vec{a}\vec{c}+\vec{d}(\vec{a}+\vec{c})=c^2+\vec{a}\vec{c}-\vec{d}(\vec{b}+\vec{d})=c^2-d^2+\vec{a}\vec{c}-\vec{b}\vec{d}=$$ $$=c^2-d^2+\frac{1}{2}(b^2+d^2-a^2-c^2)=\frac{1}{2}(b^2-d^2+c^2-a^2)$$ and we are done!