Is there a simple way to force matrix brackets to stretch according to phantom contents of `array` environment?

Here is a solution with saveboxes:

\documentclass{book}


    % ========== Table packages ==========
    \usepackage{tabularx}
% BEGIN_FOLD

        \newcommand{\arraystretchdefaultl}{1.3}
        \newcommand{\arraystretchdefault}{ \renewcommand{\arraystretch}{\arraystretchdefaultl} }
        \arraystretchdefault

% END_FOLD

    \usepackage{booktabs}
    \usepackage{makecell}


    % ========== Equations and math packages ==========
    \usepackage{amsmath}
    \usepackage{amssymb, amsfonts}
    \usepackage{mathtools}
    \usepackage[defaultmathsizes, italic, eulergreek, symbolre, symbolmisc]{mathastext}

\newsavebox{\mybox}
\newsavebox{\mysbox}
\begin{document}

\begin{align}
    \left[
        \begin{array}{c}
            \mathnormalbold{V}_{G}
            \\
        %
            \Xhline{0.25mm}
            V_{R, \thinspace i}
        %
            \\
            V_{R, \thinspace j}
        \end{array}
    \right] &=
%
    \left[
        \savebox{\mybox}{\hbox{\(\begin{array}{c !{\vrule width 0.25mm}c c}
            \bar{\mathnormalbold{Z}}_{GG} & \bar{\mathnormalbold{Z}}_{GR_{i}} & \bar{\mathnormalbold{Z}}_{GR_{j}}
        \\
            \Xhline{0.25mm}
        %
            \bar{\mathnormalbold{Z}}_{R_{i}G} & \bar{Z}_{RR, \thinspace ii} & \bar{Z}_{RR, \thinspace ij}
        %
        \\
            \bar{\mathnormalbold{Z}}_{R_{j}G} & \bar{Z}_{RR, \thinspace ij} & \bar{Z}_{R, \thinspace jj}
        \end{array}\)}}\savebox{\mysbox}{\hbox{\(\bar{Z}_{bus}\)}}
        \vphantom{\usebox{\mybox}}\hspace{\dimexpr(\dimexpr\wd\mybox/2)-\dimexpr(\wd\mysbox/2)}\usebox{\mysbox}\hspace{\dimexpr(\dimexpr\wd\mybox/2)-\dimexpr(\wd\mysbox/2)}
    \right]
    %
    \left[      
        \begin{array}{c}
            \bar{\mathnormalbold{I}}_{G}
        \\
        \Xhline{0.25mm}
        %
            \bar{I}_{R, \thinspace i}
        \\
        %
            \bar{I}_{R, \thinspace j}
        \end{array}
    \right]
\label{eq:op_imp_load_RER_int_op}
\\
\intertext{where}
\bar{\mathnormalbold{Z}}_{bus} &=
    \left[
        \begin{array}{c !{\vrule width 0.25mm}c c}
            \bar{\mathnormalbold{Z}}_{GG} & \bar{\mathnormalbold{Z}}_{GR_{i}} & \bar{\mathnormalbold{Z}}_{GR_{j}}
        \\
        \Xhline{0.25mm}
        %           
            \bar{\mathnormalbold{Z}}_{R_{i}G} & \bar{Z}_{RR, \thinspace ii} & \bar{Z}_{RR, \thinspace ij}
        \\
        %
            \bar{\mathnormalbold{Z}}_{R_{j}G} & \bar{Z}_{RR, \thinspace ij} & \bar{Z}_{R, \thinspace jj}
        \end{array}
    \right]
\label{eq:Z_bus_hat_def}
\end{align}

\end{document}

It could be better to put this functionality inside a \newcommand of course. Feel free to ask me for this if this way seems ok to you.