## 3 Proposed Method
### Objective Functions and Constraints  
An optimal balance is needed to enhance controllability and reduce RoCoF without significantly compromising response speed and damping of oscillatory modes. Therefore, two objective functions are proposed.  
The first objective, \(J_{1}\), is to maximize the controllability of the VSGs, ensuring responsive and robust operation. We use Gramian-based controllability metrics to evaluate controllability. A higher value of a suitable measure of the controllability Gramian, such as its trace, determinant, or smallest eigenvalue, implies that less energy is required to achieve the desired state transitions. The selection of Gramian-based controllability metrics for optimization purposes is shown in [21].  
The following inverse normalized functions are used to evaluate and optimize the controllability metrics:  
\[J_{\text{EV,min}}=\frac{\text{EV}_{\text{min0}}}{\text{EV}_{\text{min}}} \tag{6}\]  
\[J_{\text{EV2,min}}=\frac{\text{EV}_{\text{2 min0}}}{\text{EV}_{\text{2 min}}} \tag{7}\]  
\[J_{\text{EV,max}}=\frac{\text{EV}_{\text{max0}}}{\text{EV}_{\text{max}}} \tag{8}\]  
\[J_{\text{TR}}=\frac{\text{TR}_{0}}{\text{TR}} \tag{9}\]  
\[J_{\text{DET}}=\left(\frac{\text{DET}_{0}}{\text{DET}}\right)^{\frac{1}{2}} \tag{10}\]The expressions (6)-(10), respectively, incorporate the minimum eigenvalue of the Gramian, the sum of the two smallest eigenvalues, the maximum eigenvalue, the trace, and the determinant to the power of \(1/n\), where \(n\) is the system matrix dimension. All objective functions are normalized to the initial value of the corresponding Gramian-based metric, denoted in the subscript as 0.  
The second objective, \(J_{2}\), limits the RoCoF to stabilize the system frequency. The response of the system is analyzed through a time-domain simulation, and RoCoF is computed using the numerical five-point stencil derivative method because of its high accuracy. For a given discrete time-series data \(x(t_{k})\), the five-point stencil formula is applied as follows:  
\[\overset{\bullet}{\chi}(t_{k})\approx\frac{-x(t_{k+2})+8x(t_{k+1})-8x(t_{k-1 })+x(t_{k-2})}{12\Delta t} \tag{11}\]  
where \(\overset{\bullet}{\mathbf{x}}(t_{k})\) is the approximate first derivative of \(x(t)\) at time \(t_{k}\), \(\Delta t\) is the time step, and \(k\) is the index of the current time step.  
The determination of a system's RoCoF is given as follows. First, the time response from the initial conditions in (1) is calculated, and the time derivative of each VSG frequency is calculated using (11). Then, the maximum value of the absolute RoCoF per generator is identified. Finally, the mean of these maximum RoCoF values across all generators is used as an optimization criterion. Thus, the second objective \(J_{2}\), is given as  
\[J_{2}=\frac{\text{RoCoF}}{\text{RoCoF}_{0}}, \tag{12}\]  
where RoCoF is the actual value and RoCoF\({}_{0}\) is the initial value.  
The number of optimization parameters for the \(i\)-th VSG is two, i.e., the inertia constant \(H_{i}\) and the damping coefficient \(D_{i}\). Lower and upper bounds of all these parameters are defined so that a Pareto set is found in the following range:  
\[H_{i,min}<H_{i}<H_{i,max}\,D_{i,min}<D_{i}<D_{i,max}, \tag{13}\]  
where \(H_{i,min}\) and \(H_{i,max}\), and \(D_{i,min}\) and \(D_{i,max}\) represent the lower and upper limit for the inertia constant and damping coefficient of an \(i\)-th VSG.  
Furthermore, damping ratio constraints are applied to the oscillatory modes associated with complex conjugate eigenvalue pairs. The corresponding constraint for the oscillatory mode is given as  
\[\frac{\zeta_{0}}{\zeta}<1, \tag{14}\]  
where \(\zeta\) represents the actual damping ratio of this mode and \(\zeta_{0}\) is the initial value. These constraints ensure that the damping ratios remain greater than the initial values of the damping ratios.