Paper deep dive
Leveraging System-Level Observations to Inform Bayesian Learning of Model Parameters for Quantitative Verification
Simos Gerasimou, Xingyu Zhao
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Combining Bayesian learning and quantitative verification is a powerful toolset for analysing key quantitative properties of software systems, like reliability and response time. However, the accuracy and robustness of verification results strongly depend on the prior knowledge (PK) underlying Bayesian inference. This knowledge reflects original beliefs about the probability of events and typically depends on domain expertise. Using inaccurate or uninformative PK can negatively affect quantitative analysis, yielding incorrect verification results. Our EPIK approach tackles this important challenge by eliciting and embedding PK in quantitative verification equipped with Bayesian estimators. Unlike existing approaches that require PK on formal model transition parameters, EPIK leverages system-level properties that are directly observable and are linked to real-world semantics. EPIK formulates a twofold optimisation problem to derive the distributions of unknown transition parameters and then embeds these distributions to verify new or difficult-to-measure (elusive) properties. The detailed experimental evaluation using multiple variants of real-world case studies and diverse EPIK instantiations shows its effectiveness, flexibility and generality.
Tags
Links
- Source: https://arxiv.org/abs/2608.03489v1
- Canonical: https://arxiv.org/abs/2608.03489v1
Trouble viewing inline? Open PDF directly →
Full Text
97,457 characters extracted from source content.
Expand or collapse full text
by Leveraging System-Level Observations to Inform Bayesian Learning of Model Parameters for Quantitative Verification Simos Gerasimou Cyprus University of Technology, Cyprus simos.gerasimou@cut.ac.cy and Xingyu Zhao Wuhan University, China xingyu.zhao@whu.edu.cn (2026) Abstract. Combining Bayesian learning and quantitative verification is a powerful toolset for analysing key quantitative properties of software systems, like reliability and response time. However, the accuracy and robustness of verification results strongly depend on the prior knowledge (PK) underlying Bayesian inference. This knowledge reflects original beliefs about the probability of events and typically depends on domain expertise. Using inaccurate or uninformative PK can negatively affect quantitative analysis, yielding incorrect verification results. Our EPIK approach tackles this important challenge by eliciting and embedding PK in quantitative verification equipped with Bayesian estimators. Unlike existing approaches that require PK on formal model transition parameters, EPIK leverages system-level properties that are directly observable and are linked to real-world semantics. EPIK formulates a twofold optimisation problem to derive the distributions of unknown transition parameters and then embeds these distributions to verify new or difficult-to-measure (elusive) properties. The detailed experimental evaluation using multiple variants of real-world case studies and diverse EPIK instantiations shows its effectiveness, flexibility and generality. quantitative verification, Markov models, Bayesian learning †journalyear: 2026†copyright: c†conference: Proceedings of the 41st IEEE/ACM International Conference on Automated Software Engineering; October 12–16, 2026; Munich, Germany†booktitle: Proceedings of the 41st IEEE/ACM International Conference on Automated Software Engineering (ASE ’26), October 12–16, 2026, Munich, Germany†doi: 10.1145/3832783.3837529†isbn: 979-8-4007-2882-2/2026/10†ccs: Theory of computation Probabilistic computation 1. Introduction Modern software-controlled systems, such as robotics, environmental and healthcare monitoring applications, and cloud-based services, operate in complex and uncertain environments characterised by workloads, operational profiles, failures and resource availability that are stochastic in nature (De Lemos et al., 2013; Weyns et al., 2023). Designing and analysing the performance, dependability and other nonfunctional properties of such systems when deployed in these uncertain or, potentially, adversarial environments can be facilitated by rigorous model-based verification methods (Kwiatkowska et al., 2004). More specifically, probabilistic models, including Markov chains (Gerasimou et al., 2015; Kwiatkowska et al., 2007), queueing networks (Cortellessa et al., 2011; Balsamo et al., 2003) and stochastic Petri nets (Perez-Palacin and Merseguer, 2010; Tigane et al., 2022) can be leveraged to model and analyse the behaviour and operating environment of such systems. This paper focuses on probabilistic model checking (PMC), a formal method for automatically verifying quantitative aspects of stochastic systems (Baier and Katoen, 2008). PMC involves the construction of a Markov chain model that encodes the system’s behaviour over time, including states the system can reside in, the possible transitions between these states and information about the likelihood (probability) or timing (rate) of these transitions. Given properties of the required behaviour of these systems (e.g., reliability, performance) formally specified in temporal logic, the systematic analysis of the system model through automated algorithms enables assessing if the properties are satisfied (Baier and Katoen, 2008). Using mathematical reasoning to derive guarantees for achieving precisely defined levels of performance or efficiency is highly valuable for the design, analysis and adaptation of software-intensive cyber-physical systems (Kwiatkowska et al., 2022; Calinescu et al., 2017a, b). Recent advances improve the PMC efficiency and scalability (Filieri et al., 2012; Fang et al., 2021; Jansen et al., 2014; Evangelidis et al., 2026), enabling the analysis of more complex models and properties of software product lines (Ghezzi and Sharifloo, 2013), software architectures (Moreno et al., 2014) and cyber-physical systems (Zhao et al., 2019a). As a verification method, however, the verification results strongly depend on the veracity of the Markov model encoding the behaviour of the target system and, especially, the parameters (probability/rate) capturing the transitions between model states (Kwiatkowska et al., 2007). When these model parameters accurately reflect the current system behaviour, the verification results faithfully represent quantitative system aspects. In contrast, inaccurate or outdated model parameter values yield misleading results that unavoidably lead to incorrect engineering and adaptation decisions (Calinescu et al., 2012). Approaches to improve analysis and PMC reasoning use Bayesian methods to incorporate prior (domain) knowledge at design time (Jha et al., 2009) and update the Markov model using system observations at runtime (Epifani et al., 2009; Filieri et al., 2016; Calinescu et al., 2014). Notwithstanding their benefits, a fundamental premise underpinning these Bayesian-based approaches is that prior knowledge about the transition probability/rate between two model states can be extracted from domain experts or past similar system executions (Epifani et al., 2009). While this may be desirable, demanding domain experts to define their knowledge about model transition parameters (e.g., the probability of a fruit-picking robot transitioning between the positioning and fruit picking states) is a challenging and non-trivial problem that needs fine-grained system knowledge (Fang et al., 2022). Providing biased or uninformative values for model transition parameters yields inaccurate analysis results, leading to design decisions that can be detrimental to the system’s reliability and performance (Bishop et al., 2011). Instead of enforcing the definition of values for transition parameters of the underlying formal model, we propose employing the rich system-wide knowledge about properties of system behaviour typically possessed by domain experts (Guindon, 1990). These properties are directly observable and correspond to intuitive real-world semantics and actions. This can be expressed, for instance, as the reliability of the fruit-picking robot completing its task successfully (Fang et al., 2022) or the expected execution time of a service-based system that uses cloud-based services (Gerasimou et al., 2015). Since such system-wide knowledge embodies real-world semantics and is directly observable from past system executions, it aligns more naturally with the available domain expertise (Xie et al., 2021; Li et al., 2020). Hence, this system-wide knowledge can be analysed to extract appropriate values for the model transition parameters. Driven by this insight, we introduce EPIK, a Bayesian-based PMC approach that leverages expert knowledge of system-wide properties (termed PK-informed) for the systematic elicitation and embedding of model transition parameters. EPIK comprises the following key stages. First, during the knowledge elicitation stage, prior knowledge from domain experts, corresponding to past observations of system-wide properties (e.g., reliability, response time), each expressed as a probability distribution, initiates a multi-objective search-based elicitation problem (Coello, 2007). The synthesised Pareto set signifies model transition parameters that best approximate the probability distributions of the system-wide PK-informed properties. Then, during the knowledge embedding stage, EPIK exploits the synthesised Pareto set of transition parameters for the verification of elusive system properties, i.e., system-wide properties that are novel or rare, or those for which gathering information is risky or expensive. We evaluate EPIK on multiple variants from two software-contro-lled systems from different application domains: (1) a fruit-picking robot (FPR) (Fang et al., 2022) and (2) a service-based system from foreign exchange trading (Gerasimou et al., 2018), demonstrating the accuracy and effectiveness of EPIK in extracting and embedding prior knowledge. To the best of our knowledge, EPIK is the first tool-supported solution that considers the problem of extracting knowledge for the unknown transition parameters of a probabilistic model from system-wide properties. The real-world applicability of EPIK follows the broad adoption of similar work applying Bayesian learning in probabilistic analysis (Zhao et al., 2024; Jha et al., 2009; Epifani et al., 2009; Filieri et al., 2016), and the various models of software systems available in the repositories of probabilistic model checkers (Kwiatkowska et al., 2011; Dehnert et al., 2017). We also emphasise EPIK’s positioning with respect to downstream Bayesian learning approaches. EPIK is not a posterior inference approach that uses runtime observations. Instead, EPIK operates entirely during design-time (preparation phase) to systematically synthesise mathematically grounded prior distributions from observable, high-level system properties. The elicited distributions produced by EPIK serve as the initial priors (PKPK) for downstream Bayesian estimators (e.g., KAMI (Epifani et al., 2009)). This pipeline prevents the biased or inaccurate initialisations that frequently cause subsequent Bayesian approaches to fail or produce distorted verification results. The main contributions of our paper are: ∙ The EPIK approach to elicit and embed prior knowledge for Bayesian-based PMC that leverages PK-informed system-wide properties and enables the verification of elusive properties; ∙ An extensive EPIK evaluation on several variants of two real-world problems, for a wide variety of PK-informed properties and a set of unknown transition parameters. ∙ A prototype open-source EPIK tool and case study repository, available at https://github.com/gerasimou/EPIK. The remainder of this paper is structured as follows: Section 2 presents the required background on PMC and Bayesian learning. Section 3 introduces a running example to illustrate EPIK, which is detailed in Section 4. Section 5 describes EPIK’s implementation and evaluation. Section 6 discusses related work, and Section 7 summarises our results and suggests future research directions. 2. Preliminaries Probabilistic Model Checking. Probabilistic model checking (PMC) is a formal method for assessing quantitative properties, e.g., reliability, performance, cost, of systems exhibiting stochastic behaviour (Baier and Katoen, 2008). Markov chains capture the system’s stochastic behaviour, enabling the analysis of properties encoded as formal logic specifications. A Markov chain is a tuple M=(S,s0,δ,L)M=(S,s_0,δ,L), where S is a finite set of states, s0∈Ss_0∈ S is the initial state, δ is the state-transition function: ∙ δ:S×S→[0,1]δ:S× S→[0,1] for discrete-time Markov chains (DTMCs), with δ(si,sj)=pijδ(s_i,s_j)=p_ij giving the transition probability between states si,sj∈Ss_i,s_j∈ S, and ∑sj∈Sδ(si,sj)=1 _s_j∈ Sδ(s_i,s_j)=1; ∙ δ:S×S→ℝ≥0δ:S× S _≥ 0 for continuous-time Markov chains (CTMCs); δ(si,sj)=rijδ(s_i,s_j)=r_ij gives the rate of transition between states si,sj∈Ss_i,s_j∈ S. and L:S→2APL:S→ 2^AP is a labelling function assigning to each state a set of atomic propositions from APAP. PMC supports both DTMCs and CTMCs as modelling formalisms. DTMCs are widely used to represent systems with discrete, sequential behaviour, where transitions occur at discrete time steps. CTMCs model systems with continuous, time-dependent behaviour, where transitions can occur at any point in time. PMC can analyse properties in DTMCs/CTMCs by considering the probabilities of reaching safe/unsafe states and computing expected values (cost/rewards). Markov model states are labelled with atomic propositions that hold in those states. The properties to verify are expressed in temporal logic over these atomic propositions, e.g., probabilistic temporal tree logic (PCTL) (Bianco and de Alfaro, 1995; Hansson and Jonsson, 1994) for DTMCs and continuous stochastic logic (CSL) (Aziz et al., 1996) for CTMCs. State-of-the-art probabilistic model checkers (PRISM (Kwiatkowska et al., 2011), Storm (Dehnert et al., 2017)) implement efficient PMC algorithms for verifying Markov chain models and properties (Kwiatkowska et al., 2007). Parametric Model Checking. PMC based on DTMCs/CTMCs assumes that transition probabilities/rates are known constants, or can be estimated from existing data and experts at design time. Parametric model checking (ParaMC) (Daws, 2005) can analyse DTMCs/CTMCs whose transition probabilities/rates are specified as functions over a set of parameters, termed parametric DTMCs/CTMCs. The property under verification is given as a closed-form rational function of these parameters (Daws, 2005). This formulation brings a practical advantage by dividing the verification process into two steps. The computationally intensive symbolic analysis can be performed offline without strict computational constraints (Jansen et al., 2014). Then, the closed-form symbolic expressions can be reused later, e.g., when new knowledge or data of the parameters is learned, without rerunning the PMC every time. EPIK leverages the closed-form expressions to efficiently derive PK for the unknown model transition parameters. Bayesian Learning in Probabilistic Model Checking. Given the DTMC state sis_i, the transition to the next state follows a categorical distribution (Epifani et al., 2009). Due to the Markov property, the choice of the next state only depends on the current state, and the categorical distributions per state are independent. The categorical process of repeated transitions from state i, follows a multinomial distribution. Thus, learning the transition parameters reduces to the localised learning of k independent multinomial distributions, where k is the number of DTMC states (Gelman et al., 2014). From a Bayesian inference perspective, the posterior estimation needs a statistical model (the likelihood function) and a prior distribution (Berger et al., 1994). For the likelihood function, if we observe nijn_ij transitions from state sis_i to state sjs_j out of nin_i outgoing transitions from sis_i (termed “data” in the equations below), the binomial likelihood is (omitting the combinatorial factor that cancels out in the Bayes formula): (1) (data∣pij)=pijnij(1−pij)ni−nijPr( data p_ij)=p_ij^n_ij(1-p_ij)^n_i-n_ij Typically, for mathematical convenience in Bayesian inference, a conjugate Beta prior distribution111A Dirichlet distribution, which is a multivariate generalisation of the Beta distribution, can be used to model the prior of the i-th transition matrix row. is used for the likelihood (Epifani et al., 2009). Specifically, the canonical parameterisation222The shape parameters α\!, β are replaced by ni(0)=α+βn_i^(0)\!=\!α+β, pij(0)=α/(α+β)p_ij^(0)=α/(α+β) (2) (ni(0),pij(0))Beta (n_i^(0),p_ij^(0) ) of the Beta distribution allows an intuitive interpretation of pij(0)p_ij^(0) as the “best prior probability guess” and of ni(0)n_i^(0) as the sample size (strength) on which the prior estimation pij(0)p_ij^(0) is based on (Epifani et al., 2009; Walter and Augustin, 2009). After applying the Bayes rule and leveraging the conjugacy and canonical reparameterisation (Gelman et al., 2014), the posteriors are again a (ni(ni),pij(ni))Beta (n_i^(n_i),p_ij^(n_i) ) distribution with the updated parameters: (3) ni(ni)=ni(0)+ni,pij(ni)=ni(0)ni(0)+ni⋅pij(0)+nini(0)+ni⋅nijni n^(n_i)_i\!=\!n_i^(0)\!+\!n_i,\; p_ij^(n_i)\!=\! n_i^(0)n_i^(0)+n_i\!·\!p_ij^(0)+ n_in_i^(0)+n_i\!·\! n_ijn_i where the superscript ‘(0)(0)’ is the prior knowledge (before any observations), and the superscript ‘(ni)(n_i)’ gives the posterior parameters after observing nin_i outgoing transitions from state sis_i. Eq. (3) demonstrates that, upon observing nijn_ij transitions out of nin_i total transitions, the posterior pij(ni)p_ij^(n_i) is the weighted sum of the prior estimate pij(0)p_ij^(0) and the new observations nijni n_ijn_i (encoding the frequency of transitions from sis_i to sjs_j). The weights are proportional to ni(0)n_i^(0) (the prior sample size – strength) and nin_i (the “actual count” of data sample size). Smaller ni(0)n_i^(0) values indicate lower confidence in the priors, allowing the runtime data to influence the posteriors more. When ni(0)≃0n_i^(0)\! \!0, Eq. (3) simplifies to the Maximum Likelihood Estimation (Epifani et al., 2009). For CTMCs, where the Gamma-Poisson setup typically applies (Bernardo and Smith, 1994), Bayesian estimators for CTMC transition rates entail replacing the Beta priors with Gamma priors (Filieri et al., 2012). 3. Motivating Example We illustrate EPIK using an autonomous robot executing a fruit-picking task adapted from (Fang et al., 2022). The robot should perform the following operations: (i) position itself to collect the next fruit in its vicinity; (i) use its arm to collect the fruit; and (i) when fruit collection is unsuccessful, decide whether to retry the collection by repositioning itself or abandon the task and move to the next fruit. Fig. 1 shows the DTMC model of this robotic fruit-picking task (FPR). Starting from state s0s_0, the robot takes the appropriate position (e.g., using its onboard camera and Lidar perception devices) close to the next fruit. The positioning operation succeeds with probability p0,1p_0,1, and the robot moves to state s1s_1, where it performs the picking operation. Conversely, the positioning operation fails with probability 1−p0,11-p_0,1, the robot moves to state s3s_3, and the task execution ends (state s5s_5). From state s1s_1, the picking operation is successful with probability 1−p1,21\!-\!p_1,2, the robot moves to state s4s_4 and then terminates its execution by moving to state s5s_5. In contrast, the picking is unsuccessful with probability p1,2p_1,2, and the robot moves to state s2s_2. In this state, the robot needs to decide whether to reposition itself (state s0s_0) and retry the entire process or to abandon the picking operation (state s3s_3) and end the process (state s5s_5); this decision is shown probabilistically with transitions p2,0p_2,0 and 1−p2,01-p_2,0, respectively. The DTMC model is annotated with two reward functions, illustrated as rectangular boxes linked to states s0s_0, s1s_1 and s2s_2 (Fig. 1). The “time” reward function associates mean operation execution times t0t_0, t1t_1 and t2t_2 with the three operations performed by the robot. Likewise, an “energy” reward function associates mean energy consumption e0e_0, e1e_1 and e2e_2 with the same operations. Finally, the robot must validate the three system-level requirements from Table 1. DTMC of a robot executing a fruit-picking task; ($p_i,j$: unknown transition probabilities between states.) Figure 1. DTMC of a robot executing a fruit-picking task; (pi,jp_i,j: unknown transition probabilities between states.) Table 1. System-level requirements and formalised temporal logic properties for the fruit-picking robot. ID Description PCTL R1 What is the probability for completing the fruit-picking process successfully? =?[“success”] P_=?[F\; ``success"] R2 What is the expected time for completing the fruit-picking process? =?‘‘”[“done”] R^``time"_=?\![F ``done"] R3 What is the expected energy consumption for completing the fruit-picking process? =?‘‘”[“done”] R^``energy"_=?\![F ``done"] 4. EPIK High-level EPIK workflow showing its main stages for algebraic formulae extraction of system-level properties via parametric model checking, and knowledge elicitation and embedding via solving a multi-objective optimisation problem. Figure 2. High-level EPIK workflow showing its main stages for algebraic formulae extraction of system-level properties via parametric model checking, and knowledge elicitation and embedding via solving a multi-objective optimisation problem. 4.1. Problem Formulation Although the canonical parameterisation of priors (e.g., the Beta in Eq. (2)) provides an intuitive mechanism for incorporating prior knowledge from domain experts and historical data, it still targets individual parameters of a Markov model. Such formal abstraction makes it challenging for experts to specify and accurately express the knowledge from historical data. For instance, specifying the probability p1,2p_1,2 of failing to complete the picking operation entails a detailed understanding and expertise of the Markov model underpinning the FPR system. To address this issue, EPIK incorporates a novel concept that simplifies the elicitation of prior knowledge. Instead of requesting experts to express knowledge about model transition parameters, we ask them to express their knowledge about the observable system properties. Hence, EPIK seamlessly integrates with existing Bayesian estimators (Epifani et al., 2010; Zhao et al., 2020a; Filieri et al., 2015) and strengthens the practical use of Bayesian learning for probabilistic model checking. Definition 0 (PK-informed and Elusive Properties). We categorise system-level properties into two types: ∙ PK-informed properties for which domain experts possess (substantial) prior knowledge (PK) often gained via empirical observations of their value distributions in past executions of the same system or other systems exhibiting similar behaviour; ∙ Elusive properties signifying properties that are unusual or novel or even properties that are expensive or risky to observe, making them challenging to gather knowledge about. Domain experts may feel more comfortable and be more willing to provide prior knowledge for properties of the former type. However, acquiring an understanding of elusive properties is crucial for establishing the overall trustworthiness of the subject system. Example 0 (PK-informed versus Elusive Properties). Consider the fruit-picking robot (Section 3) and its formalised requirements (Table 1). Assume that domain experts have accumulated rich prior knowledge about requirements R1 and R2 from historical robotic missions, i.e., “based on the past (previous and similar) 200 FPR missions, we estimate that the average probability of a successful mission is 0.8 and the expected mission time is 4.75 seconds”. However, no battery usage information was recorded, signifying that no knowledge about requirement R3 is available. In this scenario, R1 and R2 are PK-informed properties, and R3 is the elusive property. We introduce the notations: X=x1,…,xKX\!\!=\!\!\x_1,...,x_K\ is the set of unknown model transition parameters, Y=y1,…,yIY\!=\!\y_1,...,y_I\ is the set of PK-informed properties and Z=z1,…,zJZ\!=\!\z_1,...,z_J\ is the set of elusive properties, with cardinality K, I and J, respectively. For each PK-informed property yi∈Yy_i\!∈\!Y, we assume that domain experts define a probability distribution PKi(yi)PK_i(y_i) as their prior knowledge. Each unknown parameter xkx_k can be estimated as a probability distribution dk(xk;θk)d_k(x_k; _k) where θk _k is the vector of parameters characterising distribution dkd_k. For instance, for a Beta distribution, θk _k is a vector comprising the two parameters from Eq. (2), used in (Epifani et al., 2009) for transition probabilities, and, similarly, a Gamma distribution for transition rates (Filieri et al., 2012). The boldface =[θ1,…,θK] θ\!=\![ _1,\!…\!, _K] defines a vector. Example 0 (PK-Informed Properties). We build on Example 4.2 where R1 and R2 are PK-informed and R3 is the elusive property, denoted as y1y_1, y2y_2 and z1z_1 respectively (i.e., Y=y1,y2Y=\y_1,y_2\ and Z=z1Z=\z_1\). Then, using the experts’ knowledge “based on previous and similar 200 fruit-picking missions, we know the average probability of a successful mission is 0.8 and the expected mission time is 4.75 seconds”, we can formalise the PK-informed properties PKR1PK_R1 and PKR2PK_R2, using the canonical parameterisation, as y1∼(200,0.8)y_1 (200,0.8) and y2∼(200,4.75)y_2 (200,4.75). Since y1y_1 is a probability and y2y_2 is a reward, we use Beta and Gamma, respectively, to match their support [0,1][0,1] and [0,+∞][0,+∞]. 4.2. Algebraic Formulae Extraction EPIK, whose high-level workflow is shown in Fig. 2, comprises three stages: first, extracting algebraic expressions for the PK-informed Y and elusive Z properties; then, eliciting knowledge about transition parameters X through the PK-informed properties Y; and, finally, embedding the derived knowledge to verify the elusive properties Z. The first EPIK stage leverages ParaMC (Section 2) to derive closed-form rational functions fif_i and gjg_j for the PK-informed Y and elusive properties Z such that yi=fi(x1,…,xK),∀i=1..Iy_i=f_i(x_1,\!…\!,x_K),∀ i=1..I, and zj=gj(x1,…,xK),∀j=1..Jz_j=g_j(x_1,\!…\!,x_K),∀ j\!=1..J. Extracting these algebraic formulae facilitates the accelerated execution of the other EPIK stages. Since we can transform the knowledge elicitation and embedding stages into optimisation problems that can be solved independently, we avoid the repeated invocation of the PMC in the loop during the optimisation process. By avoiding the expensive PMC invocation, EPIK avoids both the concrete model construction and model verification against the selected set of Y PK-informed properties for each parameter instantiation. If, however, extracting the algebraic formulae is impossible due to the complexity of the Markov model or the property (Fang et al., 2023; Jansen et al., 2014), invoking the probabilistic model checker (e.g., PRISM) iteratively would yield the same outputs, albeit much slower. Accordingly, EPIK’s algebraic formulae extraction step is recommended, as it enables faster execution of knowledge elicitation described next, but it is not mandatory. Nevertheless, if a probabilistic model checker is used in the loop, EPIK would operate equally well and yield the same output. Example 0 (FPR Formulae Extraction). Assume p1,2p_1,2 and p2,0p_2,0 of the FPR DTMC (Fig. 1) are the unknown transition parameters denoted as x1x_1 and x2x_2 (i.e. X=x1,x2X\!=\!\x_1,x_2\), with known parameters p0,1=0.9p_0,1\!=\!0.9, e0=e1=e2=1e_0\!\!=\!\!e_1\!\!=\!\!e_2\!=\!1, t0=1t_0\!\!=\!\!1, t1=2t_1\!\!=\!\!2 and t2=3t_2\!\!=\!\!3. Then, the derived algebraic formulae are: y1=fR1(x1,x2)=(0.9−0.9x1)/(1−0.9x1x2)y_1=f_R1(x_1,x_2)=(0.9-0.9x_1)/(1-0.9x_1x_2) y2=fR2(x1,x2)=(2.7x1+2.8)(1−0.9x1x2)y_2=f_R2(x_1,x_2)=(2.7x_1+2.8)(1-0.9x_1x_2) z1=gR3(x1,x2)=(1.9+0.9x1)(1−0.9x1x2)z_1=g_R3(x_1,x_2)=(1.9+0.9x_1)(1-0.9x_1x_2) The unknown transition parameters x1x_1 and x2x_2 conform to Beta distributions as commonly used (Epifani et al., 2009) with optimisable parameters =[θ1,θ2] θ=[ _1, _2] such that θ1 _1 is a vector of the two Beta distribution parameters, i.e., θ1=[n1(0),p1,2(0)] _1=[n_1^(0),p_1,2^(0)] and θ2=[n2(0),p2,0(0)] _2=[n_2^(0),p_2,0^(0)] (cf. Eq. (2)). 4.3. Knowledge Elicitation During this stage, EPIK aims to find the values θ1,..,θK _1,.., _K that best approximate the distribution of the unknown model parameters x1,…,xKx_1,…,x_K such that the distance between the properties distribution provided by domain experts and the distribution produced using θ1,..,θK _1,.., _K is minimised. Formally, (4) minθ1…θKDKL(PKi(Yi)||ℙX1∼d(x1;θ1)…XK∼d(xK;θK)fi(X1,..,XK)=Yi),∀i=1..I _ subarrayc _1\\ …\\ _K subarray\!\! \D_K\!L\! (PK_i(Y_i)||P_\!\!\! subarraycX_1 d(x_1; _1)\\ …\\ X\!_K d(x_K; _K) subarray\!\!\!f_i(X_1,..,X_K)\!\!=\!\!Y_i ),\!\!∀ i=1..I \ where DKL(P||Q)\!D_KL(P||Q)\! is the Kullback–Leibler (KL) divergence (Kullback and Leibler, 1951) of the probability distributions P and Q, often used in information theory and interpreted as the amount of information lost when using Q to approximate P. Here we present the case where the K unknown transition parameters x1,..,xKx_1,\!..\!,x_K are independent, assuming they are outgoing transition parameters from K different states of a Markov model. If, however, x1x_1 and x2x_2 are outgoing transition parameters from the same DTMC state, then they cannot be assumed to be independent. Instead a joint prior distribution like a ℎ(x1,x2,1−x1−x2;α1,α2,α3)Dirichlet(x_1,\!x_2,\!1\!\!-\!\!x_1\!\!-\!\!x_2; _1, _2, _3) (where θ1 _1 and θ2 _2 collectively become α1,α2,α3 _1, _2, _3) should be used. EPIK naturally supports such cases of optimising joint prior distributions. While KL divergence aligns naturally with information-theoretic objectives, alternative statistical divergence measures such as Jensen-Shannon (JS) divergence or Wasserstein distance can be employed to further support or refine knowledge elicitation (Cai and Lim, 2022). For instance, using the JS divergence offers a symmetric, bounded alternative that eliminates numerical instabilities in regions where the support of the transition distributions does not perfectly overlap. Alternatively, the Wasserstein distance (Earth Mover’s Distance) could be integrated to exploit the underlying geometric structure of the state space, providing stable, non-vanishing gradients even when comparing disjoint probability distributions. Performing a robustness analysis of these alternative measures with respect to the accurate extraction of the K transition parameters X, particularly under varying constraints of data sparsity or noise, is an interesting future work direction. For multiple PK-informed properties (I>1I>1), Eq. (4) produces a multi-objective optimisation problem (MOOP) that involves finding the Pareto-optimal set of solutions PSs1=∈ΘK|=[θ1,…,θK]PS_s1\!\!=\!\!\ θ∈ ^K| θ=[ _1,…, _K]\ that are on the Pareto front of the MOOP. The Pareto front is a well-established concept (Pareto, 1964) describing an optimal solution set where each solution cannot be improved further in one objective without worsening its performance in another objective. Example 0 (Elicitation Stage). Following from the previous examples, the MOOP of EPIK’s elicitation stage is given by: minn1(0),p1,2(0),n2(0),p2,0(0) _n_1^(0),p_1,2^(0),n_2^(0),p_2,0^(0) \ DKL(PKR1(Y1)||ℙx1∼(n1(0),p1,2(0))x2∼(n2(0),p2,0(0))(fR1(x1,x2)=Y1)), D_KL (PK_R1(Y_1)||P_ subarraycx_1 (n_1^(0),p_1,2^(0))\\ x_2 (n_2^(0),p_2,0^(0)) subarray(f_R1(x_1,x_2)=Y_1) ), (5) DKL(PKR2(Y2)||ℙx1∼(n1(0),p1,2(0))x2∼(n2(0),p2,0(0))(fR2(x1,x2)=Y2)) D_KL (PK_R2(Y_2)||P_ subarraycx_1 (n_1^(0),p_1,2^(0))\\ x_2 (n_2^(0),p_2,0^(0)) subarray(f_R2(x_1,x_2)=Y_2) ) \ where PKR1(Y1)PK_R1(Y_1) and PKR2(Y2)PK_R2(Y_2) are the distributions of the PK-informed properties. 1 2 Input: PKi(yi)PK_i(y_i),∀i=1…I∀ i\!=\!1\!…\!I: PK-informed properties prior knowledge fif_i, ∀i=1…I∀ i\!=\!1\!…\!I: PK-informed properties algebraic expressions dk(xk;θk)d_k(x_k;\! _k),∀k=1…K∀ k\!=\!1\!…\!K: distributions for model parameters n: sample size Output: PFs1=objectives()∈ℝI∀∈PSs1PF_s1=\objectives( θ) ^I\;∀\, θ∈ PS_s1\ PSs1=∈ΘK|=[θ1,..,θK] is Pareto-optimalPS_s1=\ θ∈ ^K| θ=[ _1,.., _K] is Pareto-optimal\ 3 4 Function ElicitKnowledge(…): 5 PFs1←∅PF_s1← , PSs1←∅PS_s1← 6 while ¬Terminate(PSs1,) Terminate(PS_s1, θ) do 7 ′←GetCandidateParamVectors(,PSs1) ← GetCandidateParamVectors( θ,PS_s1) 8 foreach ′∈′ θ ∈ do 9 OBJSθ′←∅OBJS_θ ← 10 for i=1i=1 to I do 11 sample←i∅_i← 12 for k=1k=1 to K do 13 sample←xk_x_k← Random(dk(xk;θk),n)(d_k(x_k; _k),n) 14 sample←i_i← sample⌢i_i samplexk_x_k 15 sample←yfi(samplei1,…,sampleik)_y← f_i(sample_i1,…,sample_ik) 16 fittedPKi←FitDistribution(sampleyfitted_PK_i← FitDistribution(sample_y) 17 obj←GetKL(PKi(yi),fittedPKi)obj← GetKL(PK_i(y_i),fitted_PK_i) 18 OBJSθ′←OBJSθ′⌢objOBJS_θ ← OBJS_θ obj 19 domOBJSθ′←GetDominated(OBJSθ′,PFs1)\!\!dom_OBJS_θ \!\!←\! GetDominated(OBJS_θ ,\!PF_s1) 20 if |domOBJSθ′|>0|dom_OBJS_θ |>0 then 21 PFs1=PFs1∖domOBJSθ′PF_s1\!=\!PF_s1 _OBJS_θ 22 PSs1=PSs1∖|OBJSθ∈domOBJSθ′PS_s1\!=\!PS_s1 \ θ\;|\;OBJS_θ _OBJS_θ \ 23 24 if |GetDominated(PFs1,OBJSθ′)|=0| GetDominated(PF_s1,OBJS_θ )|=0 then 25 PFs1=PFs1∪OBJSθ′PF_s1=PF_s1∪\OBJS_θ \ 26 PSs1=PSs1∪′PS_s1=PS_s1∪\ θ \ 27 28 29 PFs1,PSs1=Diversify(PFs1,PS1)PF_s1,PS_s1= Diversify(PF_s1,PS_1) 30 return PFs1,PSs1PF_s1,PS_s1 31 32 Algorithm 1 EPIK Elicitation Stage We solve the MOOP from Eq. (4) by transforming it into a search-based optimisation problem and leveraging the capabilities of evolutionary algorithms. At the core of the solution lies a multi-objective genetic algorithm driving the synthesis of the Pareto Front PFs1PF_s1 and its corresponding Pareto set PSs1PS_s1. Algorithm 1 shows the high-level steps of the solution which receives as inputs the distributions of the PK-informed system properties based on the prior knowledge from domain experts, the algebraic expressions for these properties and information related to the distribution type of each unknown model transition parameter θk _k. Then, the evolutionary loops are executed (lines 4–25) that involve the generation of candidate solutions (line 5) and the evaluation of each solution (lines 6–23) during which the solution values are used to instantiate the distribution of unknown model parameters (lines 10–13); distribution fitting occurs for each fif_i outcome (line 14) and the KL divergence is calculated for all PK-informed properties (lines 15–16). Next, the normal Pareto dominance is executed to remove inferior solutions and add non-dominated solutions (lines 17–23). The Diversify function (line 24) employs evolutionary algorithm strategies for diversity preservation to reduce the possibility of premature convergence and also to select the solutions that will participate in the next generation. Once the loop terminates, the Pareto-optimal solutions set PSs1PS_s1 and its corresponding Pareto front PFs1PF_s1 are returned. 4.4. Knowledge Embedding Given the Pareto-optimal set of solutions PSs1PS_s1 as the outcome of the elicitation stage, EPIK executes the knowledge embedding stage to establish the verification results for the elusive properties Z and select the most appropriate solution from PSs1PS_s1. To achieve this, EPIK solves another optimisation problem. If J=1J\!=\!1 (one elusive property), the problem is single-objective; otherwise (if J>1J\!>\!1), the problem is multi-objective, meaning that a reduced Pareto set PSs2⊆PSs1PS_s2 PS_s1 will be derived. Deciding the estimates (e.g., quantiles, moments) to use and solve the optimisation problem should align with the expectations of decision-makers. Similarly, the adoption of a conservative or optimistic stance prescribes whether optimising the elusive properties Z involves their maximisation or minimisation. Without loss of generality, the following formulation assumes the minimisation of the expected value for all elusive properties: (6) min[θ1,..,θK]∈PSs1X1∼d(x1;θ1)…XK∼d(xK;θK)[gj(X1,…,XK)],∀j=1..J _[ _1,.., _K]∈ PS_s1 \E_ subarraycX_1 d(x_1; _1)\\ …\\ X_K d(x_K; _K) subarray[g_j(X_1,…,X_K)],∀ j=1..J \ The algorithm for extracting the Pareto set PSs2PS_s2 and its corresponding Pareto front PFs2PF_s2 is similar to Algorithm 1. The key difference pertains to the main loop (lines 4–5) where the GetCandidateParamVectors function retrieves solutions from PSs1PS_s1 (instead of synthesising new) and the Terminate function holds when all solutions in PSs1PS_s1 have been examined. We omit this algorithm for brevity reasons. Example 0 (Embedding Stage). Given the Pareto-optimal set from the elicitation stage (Example 4.5) and elusive property R3 that will be maximised (for conservative purposes), the single-objective problem during EPIK’s embedding stage is given by: (7) max[n1(0),p1,2(0),n2(0),p2,0(0)]∈PSs1x1∼(n1(0),p1,2(0))x2∼(n2(0),p2,0(0))[gR3(x1,x2)] _[n_1^(0),p_1,2^(0),n_2^(0),p_2,0^(0)]∈ PS_s1E_ subarraycx_1 (n_1^(0),p_1,2^(0))\\ x_2 (n_2^(0),p_2,0^(0)) subarray[g_R3(x_1,x_2)] Complexity Analysis. In EPIK’s elicitation stage (Algorithm 1), since the time complexity of simple random sampling is generally O(n)O(n), lines 10–12 together yield time complexity O(n∗K)O(n*K). Since fif_i is a closed-form rational function with O(1)O(1) and each sample has a size of n, the time complexity of line 13 is O(n)O(n). Regarding FitDistribution (line 14), while the time complexity for distribution fitting can vary (it depends on the specific algorithm and the characteristics of the distribution being fitted), it is often333Specifically, in Section 5, we first calculate the mean (in O(n)O(n) time complexity) and variance (in O(n)O(n) time by calculating the sum of squared differences, in addition to the mean) of the sample with size n, then derive the two distribution parameters for Beta/Gamma from them (in O(1)O(1) constant time). Thus, the total time complexity is O(n)+O(n)+O(1)=O(n)O(n)+O(n)+O(1)=O(n). represented as O(n)O(n) where n is the number of data points used for the fitting. The time complexity of the KL divergence calculation (line 15) DKL(P||Q)=∑iP(i)logP(i)Q(i)D_KL(P||Q)= _iP(i) P(i)Q(i) is O(n)O(n) by definition, due to the need to iterate through all n elements to compute P(i)P(i) and Q(i)Q(i). Since this is executed for all I PK-informed properties, the total time complexity is for evaluating a candidate solution θ′θ is O(n∗K∗I)O(n*K*I). Concerning the evolution loop itself, this is computationally intensive and the complexity really depends on the chosen evolutionary algorithm. For example, EPIK instrumented with a multi-objective genetic algorithm has a general complexity O(I∗G∗S2)O(I*G*S^2) where I is the number of objective functions, S is the population size and G is the number of generations444A genetic algorithm like NSGA-I (Deb et al., 2002) has time complexity O(I∗S2)O(I*S^2) per generation (Sastry et al., 2005, Chapter 10.4). Over G generations, the total time complexity is O(G∗I∗S2)O(G*I*S^2).. Consequently, the time complexity of the whole Algorithm 1 to the size of various inputs is O(I2∗K∗n∗G∗S2)O(I^2*K*n*G*S^2). The key differences to note for EPIK’s knowledge embedding algorithm are: (1) the set of objectives becomes a set of statistics (that are of practical interest in the given application) on the distributions of the J elusive properties as encoded in Eq. (6); and (2) the search space of the multi-objective optimisation problem is the result of the elicitation stage, i.e., the variables to be optimised are constrained in PSs1PS_s1. Accordingly, the time complexity is O(J2∗K∗n∗G∗S2)O(J^2*K*n*G*S^2); since similar time analysis steps apply for this algorithm too, we omit the details for brevity. 5. Evaluation 5.1. Research Questions RQ1 (Accuracy): How closely can EPIK approximate the true distribution of PK-informed properties? We analyse if EPIK can derive the distributions of unknown transition parameters, thus yielding distributions of PK-informed properties that closely approximate the ground truth of those properties. RQ2 (Effectiveness): How do different optimisation approaches affect the effectiveness of EPIK? We examine the impact of different multi-objective optimisation algorithms in EPIK’s performance. We study the quality of EPIK-derived model transition parameters using the established evolutionary algorithms NSGA-I (Deb et al., 2002), SPEA2 (Zitzler et al., 2001) and CMA-ES (Igel et al., 2007). RQ3 (Conformance): How does using knowledge of PK-informed properties of varying levels of conformance affect EPIK? Deciding unknown transition parameters depends on the prior knowledge encapsulated in the set of PK-informed properties. Since the encoded prior knowledge about these properties can be conflicting, with this research question, we examine its impact on EPIK’s convergence. RQ4 (Knowledge Embedding): Can EPIK support verifying elusive properties for which prior knowledge is absent? Supporting decision-making entails providing useful insights into the trade-offs between transition parameter distributions elicited using EPIK. Thus, we explore how EPIK supports the verification of elusive properties. 5.2. Evaluation Methodology Software Systems. We evaluated EPIK using several variants of two software-controlled systems from different application domains: (1) the fruit-picking robot (FPR) (Fang et al., 2022) (Section 3); and (2) a service-based system for foreign exchange trading (FX) presented in (Gerasimou et al., 2015). FX Description (FX). An FX trader can use the system in two modes. Using the expert mode, FX iteratively analyses market activity, identifies patterns that satisfy the trader’s objectives, and automatically carries out trades. To this end, a Market watch operation extracts real-time exchange rates (bid/ask price) of selected currency pairs, which is used by a Technical analysis operation to assess the current trading conditions, predict future price movement, and decide if the trader’s objectives are: (i) “met” (causing the invocation of an Order service to carry out a trade); (i) “not met” (resulting in a new Market watch invocation); or (i) an error occurred (triggering an Alarm operation to notify the trader about discrepancies/opportunities not covered by the trading objectives). Using the normal mode, FX assesses the economic prospects of a country via a Fundamental analysis operation that collects, analyses and evaluates information (e.g., news reports, economic data and political events), and provides an assessment on the country’s outlook. If satisfied with this assessment, the trader can use the Order operation to sell/buy currency; then, a Notification operation confirms the trade completion. Finally, FX engineers are interested in analysing the system-level requirements in Table 2. Table 2. FX system-level requirements and properties ID Description PCTL R1 What is the probability that FX will complete the workflow successfully? =?[“success”] P_=?[F\; ``success"] R2 What is the expected response time for completing workflow execution? =?‘‘”[“done”] R^``time"_=?\![F ``done"] R3 What is the total cost of the third-party services used for a workflow execution? =?‘‘”[“done”] R^``cost"_=?\![F ``done"] Table 3. FPR and FX system variants analysed using EPIK Variant Details Size FPRc Y=R1; Z=R2,R3; X=p2,0p_2,0 1E+08 FPR1 Y=R1,R2; Z=R3; X=p2,0p_2,0 1E+08 FPR2 Y=R1,R2; Z=R3; X=p1,2,p2,0p_1,2,p_2,0 1E+16 FPR3 Y=R1,R2; Z=R3; X=p0,1,p1,2,p2,0p_0,1,p_1,2,p_2,0 1E+24 FX1 Y=R1,R2; Z=R3; X=pOr1p_Or1 1E+8 FX2 Y=R1,R2; Z=R3; X=pOr1,pFA2p_Or1,p_F\!A2 1E+16 FX3 Y=R1,R2; Z=R3; X=pOr1,pFA2,,pTA1p_Or1,p_F\!A2,,p_T\!A1 1E+24 Experimental Setup. We executed a diverse set of experiments using the FPR and FX system variants from Table 3. The ‘Details’ column indicates the PK-informed properties set (Y), the elusive properties set (Z) and the unknown model transition parameters (X). The ’Size’ column reports the analysis cost an exhaustive search would incur, assuming two-decimal precision for each double-valued parameter needed to encode an unknown parameter xkx_k. We equipped the optimisation algorithm of EPIK using the established multi-objective genetic algorithms NSGA-I (Deb et al., 2002) and SPEA2 (Zitzler et al., 2001) and CMA-ES (Igel et al., 2007), a state-of-the-art derivative-free evolution strategy based on covariance matrix adaptation. Unlike NSGA-I and SPEA2, that employ crossover and mutations, CMA-ES models the search space as a multivariate normal distribution and continuously learns the correlation between decision variables (parameters) and adapts its search to steer the population toward optimal regions. All algorithms leverage the DEAP Python package (De Rainville et al., 2012). We also used the following configuration to evaluate our approach: 2,000 evaluations with an initial population of 100 individuals (i.e., 20 generations in total), and default values for single-point crossover probability pc=0.9p_c=0.9 and uniform polynomial mutation probability pm=0.8p_m=0.8. We selected these values after carrying out a set of preliminary experiments and following the standard practice in the field of search-based software engineering (Gerasimou et al., 2015; Harman et al., 2012; Arcuri and Briand, 2011). For problem variants with ≥2≥\!2 objectives (PK-informed properties), i.e., |Y|>2|Y|\!>\!2, we use standard Pareto dominance (4) to construct the Pareto front approximation (Coello, 2007). To reduce the potential impact of non-deterministic behaviour affecting EPIK’s performance (e.g., when choosing the crossover point for the genetic algorithms), we adopted the established procedure in search-based software engineering (Harman et al., 2012). Thus, we performed 30 independent runs per system variant from Table 3 and each multi-objective optimisation algorithm (Arcuri and Briand, 2011). All experiments were run on a MacBook Pro with Apple M2 Max processor and 64GB of memory. Statistical Analysis. Since the actual Pareto front for the real-world systems from our experimental evaluation is unknown and its exact computation is typically impossible, we used the standard practice (Zitzler et al., 2008). For each system variant, we extract the reference front comprising the nondominated Pareto front approximation from all the runs across all EPIK instances. We used this reference front and the Pareto front quality indicators below to quantify the ‘goodness of fit’ of Pareto front approximations produced by EPIK. We use a boxplot to show the central tendency and distribution per indicator. ∙ IHVI_HV (Hypervolume): measures the objective space volume consumed by a Pareto front approximation compared to a reference front. It shows diversity and convergence, and is strictly Pareto compliant (Zitzler et al., 2008). Larger IHVI_HV values give better Pareto front approximations. ∙ IϵI_ε (Unary additive epsilon): measures the minimum value needed by a Pareto front approximation to dominate the reference front. IϵI_ε shows convergence to the reference front and is Pareto compliant. Smaller IϵI_ε values denote better Pareto front approximations. ∙ IGDI_IGD (Inverted Generational Distance): the Euclidean distance in the objective space between the Pareto front approximation and the reference front. IGDI_IGD shows diversity and convergence to the reference front. A smaller IGDI_IGD gives better Pareto approximations. Figure 3. Distributions for the PK-informed and elusive properties for the FPRc controlled variant (1 PK-informed property) and FPR1 (2 PK–informed properties) over 30 independent runs, where ground truth R1=0.80, R2=4.584 and R3=2.755. @color@gray@stroke0 @color@gray@fill0Distributions for the PK-informed and elusive properties for the FPRc controlled variant (1 PK-informed property) and FPR1 (2 PK–informed properties) over 30 independent runs, where ground truth R1=0.80, R2=4.584 and R3=2.755. Adopting the recommended practice (Arcuri and Briand, 2011), we employed inferential statistics to compare the quality indicator values obtained by EPIK variants. Thus, we confirmed that the quality indicator values do not follow a normal distribution using the Shapiro-Wilk test. Then, we used the Mann-Whitney and Kruskal-Wallis non-parametric tests with 95% confidence level (α=0.05α=0.05) to analyse the results without making assumptions about the data distribution or the variance homogeneity. Where appropriate, we did a post-hoc analysis with pairwise comparisons between the algorithms, using the conservative Bonferroni correction pcrit=α/kp_crit=α/k (k is the number of comparisons) to control the family-wise error rate. When statistical significance exists, we use Cohen’s d to quantify the importance of the observed effect (Arcuri and Briand, 2011). Cohen’s d score summarises the difference between two groups as the number of standard deviations: d=0.2d=0.2, d=0.5d=0.5 and d=0.8d=0.8 denote a small, medium and large effect size, respectively. 5.3. Results and Discussion RQ1 (Accuracy). First, we performed a controlled experiment to establish if the resulting estimates for elusive properties are accurate. Figure 3 shows the resulting distributions for the PK-informed (grey box) and elusive (yellow box) properties for the FPRc (left) and FPR1 (right) variants over 30 independent runs. The ground truth values for the elusive properties (R2=4.584 and R3=2.755 for FPRc; R3=2.755 for FPR1) were calculated using the PK-informed properties information, solving the system of linear equations to extract the expected value for the unknown model transition parameter p2,0p_2,0 and using the value for the elusive properties. Then, we used EPIK to extract values for the Beta distribution of p2,0p_2,0 and sampled from the distribution to construct the distribution of the PK-informed and elusive properties. As shown, the derived distributions closely approximate the ground truth; the delta between the distributions’ mean and (assumed) ground truth is <0.00389<0.00389, signifying EPIK’s ability to retrieve meaningful and accurate elusive property estimates. Also, we establish whether EPIK can accurately approximate the distributions of PK-informed properties by determining realistic distributions for the unknown transition parameters of the subject system models. Fig. 4 (top) shows example Pareto fronts for the FPR variants FPR1, FPR2 and FPR3 generated during the elicitation stage of EPIK. Irrespective of the evolutionary algorithm underpinning our approach (NSGA-I or SPEA2), EPIK is capable of producing Pareto front approximations with diversified KL divergence values (DKLD_KL) for the considered PK-informed properties R1 and R2. In particular, the Pareto fronts for FPR1 comprise many non-dominated solutions, demonstrating smooth end-to-end coverage of the objective space. Thus, they yield diverse parameter value pairs for the unknown transition parameters and offer decision-makers a detailed overview of the tradeoffs between the system requirements. Pareto fronts (top) and sample distributions of requirements R1 (middle) and R2 (bottom) estimated by EPIK for different numbers of unknown model transition parameters of the fruit-picking robot use case (variants FPR1, FPR2, FPR3) Figure 4. Pareto fronts (top) and sample distributions of requirements R1 (middle) and R2 (bottom) estimated by EPIK for different numbers of unknown model transition parameters of the fruit-picking robot use case (variants FPR1, FPR2, FPR3) The Pareto front approximations produced for FPR2 and FPR3, albeit sparser, also cover the objective space satisfactorily, with the majority of the solutions clustered at the bottom left of the graph (where the total DKLD_KL is around 1). As expected, however, the range of DKLD_KL values per PK-informed property increases with solutions exceeding 3.5 (2.0) and 6 (2.5) for property R1 (R2) for FPR2 and FPR3, respectively. Increasing the number of unknown transition parameters in FPR instances FPR2 and FPR3 unavoidably expands the search space (Table 3) considerably and reshapes the landscape of the objective space. Consequently, EPIK identifies solutions at the extremes of the objective space with very small DKLD_KL for one property and very large DKLD_KL for another. We visualise the distribution of PK-informed properties R1 (Fig. 4 middle) and R2 (Fig. 4 bottom) using EPIK-derived Pareto front solutions against the prior knowledge provided by domain experts for these properties. Evidently, the EPIK-produced distributions are good approximations of the target distributions. The FPR1 solution produces distributions that accurately match the target distributions (total DKLD_KL=0.08), while the FPR2 and FPR3 are close enough but with higher total DKL≈D_KL\!\! 0.27 (given the extra ’degrees of freedom’). We obtained similar Pareto front approximations and sampled distributions for the PK-informed properties using EPIK for the FX problem instances FX1, FX2 and FX3; due to space constraints, these results are available on the project webpage. These results corroborate our findings and demonstrate EPIK’s capabilities to extract combinations of unknown transition parameter values that closely approximate the properties distributions given by domain experts. KL divergence between EPIK solutions and Random search against ground truth for FPRc (left); and quality indicator boxplots for NSGAII and Random for FPR variants (right) Figure 5. KL divergence between EPIK solutions and Random search against ground truth for FPRc (left); and quality indicator boxplots for NSGAII and Random for FPR variants (right) Boxplots comparing EPIK with multi-objective evolutionary algorithms for the FPR (top) and FX (bottom) systems Figure 6. Boxplots comparing EPIK with multi-objective evolutionary algorithms for the FPR (top) and FX (bottom) systems RQ2 (Effectiveness). First, Fig. 5 shows the KL result between EPIK and Random search against ground truth for FPRc (left); and quality indicator boxplots for NSGAII and Random for FPR variants (right). These results show a clear gap between EPIK-derived solutions and Random, signifying that the optimisation problem addressed by EPIK is non-trivial, mandating the use of sophisticated evolutionary algorithms underpinning EPIK. Next, we investigate EPIK’s effectiveness in eliciting admissible values for encoding the distributions of the unknown transition parameters using diverse multi-objective optimisation algorithms. Fig. 6 depicts boxplots, over 30 independent runs, of the quality indicators IϵI_ε, IGDI_IGD and IHVI_HV for the various FPR and FX system instances when EPIK is instrumented with the widely-used evolutionary algorithms NSGA-I (Deb et al., 2002), SPEA2 (Zitzler et al., 2001) and CMA-ES (Igel et al., 2007). We observed that, despite EPIK’s flexibility, the CMA-ES-based EPIK underperformed compared to the other EPIK instances across all problem variants and for all quality indicators. The statistical analysis using Kruskal-Wallis yielded a statistically significant difference for all problem variants and quality indicator combinations (p-value << 0.0002), substantiating the performance difference that is evident in the boxplots. Through the post-hoc analysis, involving pairwise comparisons using Mann-Whitney and Cohen’s d effect size, we established a statistically significant difference (p-value << 1.205E-10) with a very large effect size (d>>1.21) for all NSGA-I and SPEA2 comparisons against CMA-ES. In contrast, we found no statistically significant difference between EPIK using NSGA-I or SPEA2; the p-value for all problem variants and indicators was [0.051, 0.79], exceeding the significance level α=0.05α=0.05. These findings clearly evidence that EPIK employing NSGA-I and SPEA2 can produce Pareto front approximations that yield significantly better quality indicators than CMA-ES. Furthermore, these results establish EPIK’s generality in identifying effective distributions for the unknown transition parameters by leveraging different multi-objective evolutionary algorithms. Distribution of EPIK results for various conformance levels of FPR PK-informed properties Figure 7. Distribution of EPIK results for various conformance levels of FPR PK-informed properties RQ3 (Conformance). We examine EPIK’s convergence when the PK-informed properties incorporate conflicting prior knowledge. This conflict can occur when domain experts have differing opinions about the PK-informed properties, leading to divergent accumulated knowledge. A similar conflicting situation can occur when the same expert has strong knowledge about one PK-informed property but weak or contradictory knowledge about the other properties. We investigated this situation through four different conformance scenarios for the PK-informed properties of system FPR1. Each scenario encodes a different conformance level, ranging from high conformance to high conflict. Fig. 7 shows the resulting distributions of the best solution with the lowest total DKLD_KL (given equal weight to objectives) found by the NSGA-I-based EPIK. The solutions produced by EPIK for the conformance and low-conflict scenarios yielded distributions that closely approximate the distributions of PK-informed properties from domain experts. Unavoidably, the higher the conflict level, the more divergent the distributions (especially for the medium and high conflict scenarios), indicating that such conflicting knowledge prevents EPIK from finding solutions that further reduce the total DKLD_KL values. The boxplots in Fig. 8 (left), generated over 30 independent runs, clearly show the distance between the total DKLD_KL values for the different conformance scenarios. Through the execution of the Kruskal-Wallis and Mann-Whitney inferential tests, we established statistically significant differences (p≪\; \;0.05 and very large effect size) for all pairwise comparisons, including the conformance scenarios. The narrow width of all boxplots confirms that the resulting gap is not due to randomness in EPIK’s generalisation or operation. Lastly, Fig. 8 (right) shows the median DKLD_KL evolution for each population per generation for a specific EPIK run across all conformance scenarios. The evident fluctuation in the medium and high conflict scenarios illustrates the difficulty in finding stable Pareto front approximations that yield low DKLD_KL values. This observation is corroborated by the evident gap against the conformance and low conflict scenarios. Boxplots for conformance comparison on FPR (left) and KL Evolution (median per generation) for a specific run and various conformance levels (right) Figure 8. Boxplots for conformance comparison on FPR (left) and KL Evolution (median per generation) for a specific run and various conformance levels (right) These results provide strong empirical evidence that EPIK’s ability to synthesise effective distributions for the unknown parameter values depends on the conformance level of the prior knowledge of the PK-informed properties. Identifying conflicting prior knowledge would be an interesting direction for future work to avoid wasting resources and to involve humans in the loop to resolve such conflicts. Elusive property R3 distribution for 3 FPR1 Pareto-optimal solutions from EPIK’s knowledge embedding stage Figure 9. Elusive property R3 distribution for 3 FPR1 Pareto-optimal solutions from EPIK’s knowledge embedding stage RQ4 (Knowledge Embedding). We employed the Pareto-optimal approximation solutions set synthesised during EPIK’s elicitation stage to perform knowledge embedding by solving the optimisation problem in Eq. (6) and quantifying elusive property R3 (Table 1). Fig. 9 shows the distributions of the verification results of elusive property R3 and the estimated quantiles q=0.05,0.5,0.95q=\0.05,0.5,0.95\ for three FPR1 Pareto-optimal solutions p2,0=Beta(8.07,1.6)p_2,0\!\!=\!\!Beta(8.07,1.6), p2,0=Beta(9.83,1.75)p_2,0\!\!=\!\!Beta(9.83,1.75) and p2,0=Beta(10.67,3.16)p_2,0\!\!=\!\!Beta(10.67,3.16). The different Pareto-optimal solutions (signifying different p2,0p_2,0 instantiations) enable deriving elusive property distributions, which, albeit looking similar, yield different quantile values and density shapes. Decision-makers can leverage this information to examine the shape of the elusive property’s distribution. Since no information existed for the unknown transition parameter values of the considered system, without EPIK’s embedding stage, this analysis could not have been performed at all or would involve a biased guess of the transition parameters, incurring the risk of inaccurate outcomes. EPIK instead leverages the prior knowledge about the PK-informed properties to perform a well-reasoned analysis, enabling decision-makers also to select, as per Eq. (6), an appropriate solution that conforms to the system and elusive property semantics. For instance, adopting a conservative approach for quantile q=0.5q=0.5, and since the elusive property R3 quantifies the expected energy consumption for completing the fruit-picking process, decision-makers would select the middle solution p2,0=Beta(9.83,1.75)p_2,0\!\!=\!\!Beta(9.83,1.75) that yields the maximum median value. The selected solution can then drive runtime quantitative verification (Filieri and Tamburrelli, 2013; Zhao et al., 2024) and instrument self-adaptation (Calinescu et al., 2017a; Filieri et al., 2016). EPIK Discussion & Guidelines. 1) Encoding summary statistics as distributions: When experts provide only a summary statistic (e.g., “reliability is 0.8”), EPIK maps it to a canonical Beta/Gamma parameterisation (Eq. (2)). The point estimate serves as the “best guess” (pi,j(0)p_i,j^(0) or Gamma mean), while confidence is encoded via the sample size ni(0)n_i^(0). Setting ni(0)≈1n_i^(0)≈ 1 yields a near-uninformative prior, allowing subsequent Bayesian estimators (Epifani et al., 2009) to be driven primarily by runtime data; setting ni(0)n_i^(0) to historical observations (e.g., ni(0)=200n_i^(0)=200 for 200 past missions) encodes a strong, concentrated belief. This supports a continuous spectrum from uninformative to highly informative priors without requiring experts to specify a distribution shape by hand. 2) Conflict detection: Because EPIK’s elicitation stage (Section 4.3) minimises KL divergence (Eq. (4)), the Pareto front diagnoses conflicting prior knowledge. Adopters can monitor two indicators (RQ3): (i) the magnitude of total DKLD_KL achieved by Pareto-optimal solutions, where values clustering far from 0 indicate that no candidate distribution can simultaneously satisfy all properties; and (i) the generational stability of the median DKLD_KL (Fig. 8, right), where sustained oscillation rather than smooth convergence signals contradictory objectives. Either symptom prompts users to reconcile the PK-informed properties with experts before trusting the resulting priors. 3) Managing Non-Identifiability and Epistemic Ambiguity: In inverse problems, multiple distinct transition parameter distributions can yield identical system-level property distributions. EPIK treats this non-identifiability as a feature: the synthesised Pareto set (PSs1PS_s1) is not a single “correct” parameter recovery, nor is any θ∈PSs1θ∈ PS_s1 the absolute ground truth. The Pareto set represents the family of admissible priors consistent with the formal model and domain knowledge. 4) Conditional Soundness of Elusive Properties: Elusive properties cannot be validated directly at design time due to missing prior data. We stress that EPIK’s soundness in verifying these properties is strictly conditional: it assumes that PK-informed properties are accurate and the underlying model structure is correct. If expert-provided priors are biased, the resulting parameters, and consequently, the elusive property evaluations, will also be biased. This is an unavoidable epistemic limitation of all Bayesian elicitation techniques. 5) Performance & Scalability Trade-offs. Practitioners should configure EPIK based on the model’s complexity. Offline algebraic formula extraction has a high upfront cost (seconds to ∼8 8 minutes for the systems evaluated) but enables millisecond-level evaluations; this is recommended for DTMC/CTMC models where closed-form rational functions can be symbolically computed. If model scale or property complexity prevents symbolic extraction, adopters should fall back to iterative PMC, which guarantees correct results but incurs high continuous overhead by running PMC in-the-loop. 5.4. Threats to Validity We limit construct validity threats that may occur due to simplifications and assumptions in the evaluation by using formal models and properties of systems from real-world case studies from the literature (e.g., FPR (Fang et al., 2022), FX (Gerasimou et al., 2015)). Also, EPIK employs KL for comparing the candidate and ground truth distributions, thus leveraging its asymmetric property and optimising the candidate distribution to the ground truth. Experimenting with other distance metrics (e.g., Jensen–Shannon divergence) would help in further validating EPIK. We reduce internal validity threats that may introduce bias in establishing cause-effect relationships in our experiments by reporting results over 30 independent runs per experiment, and using statistical tests to check for statistical significance. Also, since EPIK depends on the prior knowledge from experts for PK-informed properties, we examine (RQ3) the EPIK’s ability to cope with inaccurate and/or conflicting prior knowledge. Finally, we enable replication by making all experimental results publicly available on our project webpage. We limit external validity threats that may reduce the EPIK’s generalisability by leveraging knowledge from experts to build the PK-informed properties. The strength of prior knowledge can be easily encoded into the PK-informed properties as well as the incorporation of prior knowledge from multiple experts (i.e., the cumulative knowledge encoded in the parameters of the corresponding distribution, e.g., Gamma) for the PK-informed property. Although we used established case studies, doing more experiments in domains and applications with different characteristics than those from our evaluation would further validate EPIK’s applicability and scalability. 6. Related Work Applications of quantitative verification in software engineering. PMC-based quantitative verification has been widely applied in various domains, including self-adaptive systems (Calinescu et al., 2012, 2017c) and, more recently, autonomous systems (Kwiatkowska et al., 2022; Vázquez et al., 2025). Accordingly, there has been a significant focus of research and publications within the software engineering community. The QoSMOS framework (Calinescu et al., 2011) is an illustrative example for developing adaptive service-based systems, using PRISM to compute quality-of-service properties. Similarly, the application in (Filieri and Tamburrelli, 2013) is a typical web-based self-adaptive system that comprises an HTTP Proxy server, a web server and an application server. A low-power wireless bus system is studied in (Filieri et al., 2016) to be efficiently verified at run time as soon as changes occur, while dynamic power management systems are also investigated (Calinescu et al., 2017a; Gerasimou et al., 2015). Interested readers can find further examples of PMC applications in software engineering for self-adaptive systems in (Vogel, ). For autonomous systems, PMC-based quantitative verification has been applied to a range of scenarios, including spacecraft reconfiguration (Nardone et al., 2016), motion planning (Lahijanian et al., 2011, 2015), controller synthesis for unmanned underwater vehicles and unmanned aerial vehicles (Giaquinta et al., 2018; Getir Yaman et al., 2025; Zhao et al., 2019a; Calinescu et al., 2015), and safety/reliability assurance in harsh environments (Calinescu et al., 2017c; Zhao et al., 2019b; Dong et al., 2022; Zhao et al., 2024; Gerasimou et al., 2021). Additional applications include task allocation and planning for mobile robots (Lacerda et al., 2019; Getir Yaman et al., 2025; Vázquez et al., 2026) and their battery charge scheduling (Zhao et al., 2019a; Tomy et al., 2020). While this is not an exhaustive list of practical applications for quantitative verification, all of these scenarios can benefit from EPIK when integrated with Bayesian learning methods to develop more accurate probabilistic models. Parameter estimation in quantitative verification. Despite progress in quantitative verification (Kwiatkowska et al., 2022), a key challenge persists: quantitative verification assumes that the models accurately reflect real software behaviour. This is often true for model structures, but transition probabilities/rates are harder to estimate correctly (Calinescu et al., 2012). Relying on point estimates provided by domain experts, inferred through model fitting (Su and Rosenblum, 2013) or updated at runtime (Filieri et al., 2011) incurs unquantified estimation errors, which are propagated and compounded in later verification steps in ways that are unknown but likely to be significant (Calinescu et al., 2012). Although recent research (Calinescu et al., 2016) synthesises bounds for unknown transition parameters via confidence intervals, it considers solely operational data, disregarding any prior human knowledge. This limitation motivates the use of Bayesian estimators (Filieri et al., 2011). Bayesian learning for runtime quantitative verification. The KAMI framework (Epifani et al., 2009) introduced Bayesian learning to estimate the transition probabilities of DTMCs. The framework was later retrofitted for CTMCs (Filieri et al., 2012) and extended with ageing factors to capture time-varying transition probabilities (Calinescu et al., 2014) and with a lightweight adaptive filter to reduce the noise (Filieri et al., 2015). Robust Bayesian estimators yielding interval estimates for DTMCs and CTMCs have been introduced in (Zhao et al., 2019b, 2024), while (Epifani et al., 2010; Zhao et al., 2020a) devised Bayesian-based change-point detectors to determine change points during the system operation. EPIK is unique as it facilitates the elicitation and embedding of prior knowledge required by these techniques. Domain knowledge elicitation and embedding. This topic is studied across software engineering (O’Hagan et al., 2006; Wright and Ayton, 1987), including requirements engineering (Nuseibeh and Easterbrook, 2000; Hadar et al., 2014), where knowledge can help clarify the needs of stakeholders and the system boundaries so as to improve software design (Adelson and Soloway, 1986; Sonnentag, 1998). Likewise, in software testing, domain knowledge aids in test case prioritisation and selection (Jin et al., 2019), while leveraging prior knowledge in safety-critical systems testing can significantly accelerate the testing process, even when the PK is partial and vague (Zhao et al., 2020b). Incorporating domain prior knowledge is also gaining attention in machine learning testing and verification (Braiek and Khomh, 2020; Stewart and Ermon, 2017; Xie et al., 2021; Kerrigan et al., 2021), backdoor attack/defence (Huang et al., 2022), and explainability (Heaton and Fung, 2023). To the best of our knowledge, only (Oghabi et al., 2011; Zervoudakis et al., 2013; Calinescu et al., 2021) combine domain knowledge with PMC-based quantitative verification. The Web Ontology Language is used to describe service behaviours as domain knowledge and then generate stochastic formal models (Oghabi et al., 2011). Similarly, a YAML-based domain-specific language is used to express system specifications that can be compiled into PCTL properties and verified by a model checker (Zervoudakis et al., 2013). Domain-specific modelling patterns in (Calinescu et al., 2021) and fragmentation techniques (Fang et al., 2023) speed up the analysis of parametric Markov chains, enabling the extraction of closed-form expressions for properties of such patterns. This research devises formal models and properties exploiting other explicit formalisms. EPIK is unique and innovative as it is natively designed to elicit and embed prior knowledge into model parameters without relying on external semantics or domain-specific languages. 7. Conclusion & Future Work We introduced EPIK, the first work for eliciting and embedding prior knowledge in Bayesian learning for quantitative verification. EPIK leverages expert knowledge of observable system-level properties to formulate a two-fold optimization problem, extracting Pareto-optimal approximation sets of unknown transition parameters whose Pareto fronts closely match expert expectations. Decision-makers can use these fronts to quantify elusive (rare, novel, or expensive) system properties. Through a comprehensive experimental evaluation using multiple variants of real-world case studies and diverse EPIK instantiations, we demonstrated the effectiveness, flexibility and generality of our approach. Our future work includes (1) enabling knowledge elicitation and embedding for interval Bayesian verification (Zhao et al., 2024); (2) investigating parallelisation methods to improve EPIK’s scalability; (3) evaluating EPIK to other case studies and scenarios. Data Availability. The open-source EPIK replication package is available at https://doi.org/10.5281/zenodo.19337648. References B. Adelson and E. Soloway (1986) A model of software design. International Journal of Intelligent Systems 1 (3), p. 195–213. Cited by: §6. A. Arcuri and L. Briand (2011) A practical guide for using statistical tests to assess randomized algorithms in software engineering. In 2011 33rd International Conference on Software Engineering (ICSE), p. 1–10. Cited by: §5.2, §5.2, §5.2. A. Aziz, K. Sanwal, V. Singhal, and R. Brayton (1996) Verifying continuous time Markov chains. In Computer Aided Verification, R. Alur and T. A. Henzinger (Eds.), LNCS, Vol. 1102, p. 269–276. External Links: ISBN 978-3-540-68599-9 Cited by: §2. C. Baier and J. Katoen (2008) Principles of model checking. MIT press. Cited by: §1, §2. S. Balsamo, V. D. N. Personè, and P. Inverardi (2003) A review on queueing network models with finite capacity queues for software architectures performance prediction. Performance Evaluation 51 (2-4), p. 269–288. Cited by: §1. J. O. Berger, E. Moreno, L. R. Pericchi, M. J. Bayarri, J. M. Bernardo, J. A. Cano, J. De la Horra, J. Martín, D. Ríos-Insúa, B. Betrò, et al. (1994) An overview of robust Bayesian analysis. Test 3 (1), p. 5–124. Cited by: §2. J. M. Bernardo and A. F. M. Smith (1994) Bayesian theory. Wiley (eng). External Links: ISBN 0-471-92416-4 Cited by: §2. A. Bianco and L. de Alfaro (1995) Model checking of probabilistic and nondeterministic systems. In Foundations of Software Technology and Theoretical Computer Science, P. S. Thiagarajan (Ed.), LNCS, Vol. 1026, p. 499–513. External Links: ISBN 978-3-540-49263-4 Cited by: §2. P. Bishop, R. Bloomfield, B. Littlewood, A. Povyakalo, and D. Wright (2011) Toward a formalism for conservative claims about the dependability of software-based systems. IEEE Transactions on Software Engineering 37 (5), p. 708–717. External Links: Document Cited by: §1. H. B. Braiek and F. Khomh (2020) On testing machine learning programs. Journal of Systems and Software 164, p. 110542. Cited by: §6. Y. Cai and L. Lim (2022) Distances between probability distributions of different dimensions. IEEE Transactions on Information Theory 68 (6), p. 4020–4031. Cited by: §4.3. R. Calinescu, L. Grunske, M. Kwiatkowska, R. Mirandola, and G. Tamburrelli (2011) Dynamic QoS management and optimization in service-based systems. IEEE Transactions on Software Engineering 37 (3), p. 387–409. External Links: ISSN 0098-5589 Cited by: §6. R. Calinescu, M. Autili, J. Cámara, A. Di Marco, S. Gerasimou, P. Inverardi, A. Perucci, N. Jansen, J. Katoen, M. Kwiatkowska, et al. (2017a) Synthesis and verification of self-aware computing systems. In Self-Aware Computing Systems, p. 337–373. Cited by: §1, §5.3, §6. R. Calinescu, S. Gerasimou, and A. Banks (2015) Self-adaptive software with decentralised control loops. In Fundamental Approaches to Software Engineering, A. Egyed and I. Schaefer (Eds.), LNCS, Vol. 9033, p. 235–251. External Links: ISBN 978-3-662-46675-9 Cited by: §6. R. Calinescu, S. Gerasimou, K. Johnson, and C. Paterson (2017b) Using runtime quantitative verification to provide assurance evidence for self-adaptive software. In Software Engineering for Self-Adaptive Systems I. Assurances, p. 223–248. Cited by: §1. R. Calinescu, C. Ghezzi, K. Johnson, P. Mauro, Y. Rafiq, and G. Tamburrelli (2016) Formal verification with confidence intervals to establish quality of service properties of software systems. IEEE Transactions on Reliability 65 (1), p. 107–125. External Links: ISSN 0018-9529 Cited by: §6. R. Calinescu, C. Ghezzi, M. Kwiatkowska, and R. Mirandola (2012) Self-adaptive software needs quantitative verification at runtime. Communication of the ACM 55 (9), p. 69–77. External Links: ISSN 0001-0782 Cited by: §1, §6, §6. R. Calinescu, C. Paterson, and K. Johnson (2021) Efficient Parametric Model Checking Using Domain Knowledge. IEEE Transactions on Software Engineering 47 (6), p. 1114–1133. External Links: Document Cited by: §6. R. Calinescu, Y. Rafiq, K. Johnson, and M. E. Bakır (2014) Adaptive model learning for continual verification of non-functional properties. In Proc. of the 5th ACM/SPEC Int. Conference on Performance Engineering, ICPE ’14, New York, NY, USA, p. 87–98. Cited by: §1, §6. R. Calinescu, D. Weyns, S. Gerasimou, M. U. Iftikhar, I. Habli, and T. Kelly (2017c) Engineering trustworthy self-adaptive software with dynamic assurance cases. IEEE Transactions on Software Engineering 44 (11), p. 1039–1069. Cited by: §6, §6. C. A. C. Coello (2007) Evolutionary algorithms for solving multi-objective problems. Springer. Cited by: §1, §5.2. V. Cortellessa, A. Di Marco, and P. Inverardi (2011) Model-based software performance analysis. Vol. 980, Springer. Cited by: §1. C. Daws (2005) Symbolic and parametric model checking of discrete-time markov chains. In Theoretical Aspects of Computing, p. 280–294. External Links: ISBN 978-3-540-31862-0 Cited by: §2. R. De Lemos, H. Giese, H. A. Müller, M. Shaw, J. Andersson, M. Litoiu, B. Schmerl, G. Tamura, N. M. Villegas, T. Vogel, et al. (2013) Software engineering for self-adaptive systems: a second research roadmap. In Software Engineering for Self-Adaptive Systems I: International Seminar, Dagstuhl Castle, Germany, October 24-29, 2010 Revised Selected and Invited Papers, p. 1–32. Cited by: §1. F. De Rainville, F. Fortin, M. Gardner, M. Parizeau, and C. Gagné (2012) Deap: a python framework for evolutionary algorithms. In Proceedings of the 14th Annual Conference Companion on Genetic and Evolutionary Computation, p. 85–92. Cited by: §5.2. K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan (2002) A fast and elitist multiobjective genetic algorithm: NSGA-I. IEEE Transactions on Evolutionary Computation 6 (2), p. 182–197. Cited by: §5.1, §5.2, §5.3, footnote 4. C. Dehnert, S. Junges, J. Katoen, and M. Volk (2017) A Storm is coming: A modern probabilistic model checker. In Computer Aided Verification, R. Majumdar and V. Kunčak (Eds.), LNCS, Vol. 10427, Cham, p. 592–600. External Links: ISBN 978-3-319-63390-9 Cited by: §1, §2. Y. Dong, X. Zhao, and X. Huang (2022) Dependability analysis of deep reinforcement learning based robotics and autonomous systems through probabilistic model checking. In 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), p. 5171–5178. Cited by: §6. I. Epifani, C. Ghezzi, R. Mirandola, and G. Tamburrelli (2009) Model evolution by run-time parameter adaptation. In the 31st Int. Conference on Software Engineering, ICSE’09, Washington, DC, USA, p. 111–121. Cited by: §1, §1, §1, §2, §2, §2, §2, §4.1, Example 4.4, §5.3, §6. I. Epifani, C. Ghezzi, and G. Tamburrelli (2010) Change-point detection for black-box services. In Proc. of the 18th ACM SIGSOFT Int. Symp. on Foundations of Software Engineering, FSE’10, New York, NY, USA, p. 227–236. Cited by: §4.1, §6. A. Evangelidis, G. Vázquez, and S. Gerasimou (2026) Accelerating policy synthesis in large-scale mdps via hierarchical adaptive refinement. Proceedings of the ACM on Software Engineering 3 (FSE), p. 4459–4481. Cited by: §1. X. Fang, R. Calinescu, S. Gerasimou, and F. Alhwikem (2021) Fast parametric model checking through model fragmentation. In 2021 IEEE/ACM 43rd Int. Conference on Software Engineering (ICSE), p. 835–846. Cited by: §1. X. Fang, R. Calinescu, S. Gerasimou, and F. Alhwikem (2023) Fast parametric model checking with applications to software performability analysis. IEEE Transactions on Software Engineering 49 (10), p. 4707–4730. Cited by: §4.2, §6. X. Fang, R. Calinescu, C. Paterson, and J. Wilson (2022) PRESTO: predicting system-level disruptions through parametric model checking. In Proceedings of the 17th Symposium on Software Engineering for Adaptive and Self-Managing Systems, p. 91–97. Cited by: §1, §1, §1, §3, §5.2, §5.4. A. Filieri, G. Tamburrelli, and C. Ghezzi (2016) Supporting self-adaptation via quantitative verification and sensitivity analysis at run time. IEEE Transactions on Software Engineering 42 (1), p. 75–99. External Links: ISSN 0098-5589 Cited by: §1, §1, §5.3, §6. A. Filieri, C. Ghezzi, and G. Tamburrelli (2011) Run-time efficient probabilistic model checking. In Proc. of the 33rd International Conference on Software Engineering, ICSE ’11, New York, NY, USA, p. 341–350. External Links: ISBN 978-1-4503-0445-0, Document Cited by: §6. A. Filieri, C. Ghezzi, and G. Tamburrelli (2012) A formal approach to adaptive software: Continuous assurance of non-functional requirements. Formal Aspects of Computing 24 (2), p. 163–186. External Links: ISSN 1433-299X Cited by: §1, §2, §4.1, §6. A. Filieri, L. Grunske, and A. Leva (2015) Lightweight adaptive filtering for efficient learning and updating of probabilistic models. In Proc. of the 37th Int. Conference on Software Engineering, ICSE ’15, Florence, Italy, p. 200–211. External Links: ISBN 978-1-4799-1934-5 Cited by: §4.1, §6. A. Filieri and G. Tamburrelli (2013) Probabilistic verification at runtime for self-adaptive systems. In Assurances for Self-Adaptive Systems: Principles, Models, and Techniques, J. Cámara, R. de Lemos, C. Ghezzi, and A. Lopes (Eds.), LNCS, Vol. 7740, p. 30–59. External Links: ISBN 978-3-642-36249-1 Cited by: §5.3, §6. A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2014) Bayesian data analysis. 3rd edition, CRC press. Cited by: §2, §2. S. Gerasimou, R. Calinescu, and G. Tamburrelli (2018) Synthesis of probabilistic models for quality-of-service software engineering. Automated Software Engineering 25 (4), p. 785–831. Cited by: §1. S. Gerasimou, J. Cámara, R. Calinescu, N. Alasmari, F. Alhwikem, and X. Fang (2021) Evolutionary-guided synthesis of verified pareto-optimal mdp policies. In 2021 36th IEEE/ACM Int. Conference on Automated Software Engineering (ASE), p. 842–853. Cited by: §6. S. Gerasimou, G. Tamburrelli, and R. Calinescu (2015) Search-based synthesis of probabilistic models for Quality-of-Service software engineering. In 2015 30th IEEE/ACM Int. Conference on Automated Software Engineering, ASE’15, Lincoln, NE, USA, p. 319–330. Cited by: §1, §1, §5.2, §5.2, §5.4, §6. S. Getir Yaman, P. Ribeiro, A. Cavalcanti, R. Calinescu, C. Paterson, and B. Townsend (2025) Specification, validation and verification of social, legal, ethical, empathetic and cultural requirements for autonomous agents. Journal of Systems and Software. Cited by: §6. C. Ghezzi and A. M. Sharifloo (2013) Model-based verification of quantitative non-functional properties for software product lines. Information and Software Technology 55 (3), p. 508–524. Cited by: §1. R. Giaquinta, R. Hoffmann, M. Ireland, A. Miller, and G. Norman (2018) Strategy synthesis for autonomous agents using prism. In 10th NASA Formal Methods Symposium (NFM’18), p. 220–236. Cited by: §6. R. Guindon (1990) Knowledge exploited by experts during software system design. International Journal of Man-Machine Studies 33 (3), p. 279–304. External Links: ISSN 0020-7373 Cited by: §1. I. Hadar, P. Soffer, and K. Kenzi (2014) The role of domain knowledge in requirements elicitation via interviews: an exploratory study. Requirements Engineering 19, p. 143–159. Cited by: §6. H. Hansson and B. Jonsson (1994) A logic for reasoning about time and reliability. Formal Aspects of Computing 6 (5), p. 512–535. External Links: ISSN 1433-299X Cited by: §2. M. Harman, S. A. Mansouri, and Y. Zhang (2012) Search-based software engineering: trends, techniques and applications. ACM Computing Surveys (CSUR) 45 (1), p. 1–61. Cited by: §5.2, §5.2. H. Heaton and S. W. Fung (2023) Explainable ai via learning to optimize. Scientific Reports 13 (1), p. 10103. Cited by: §6. W. Huang, X. Zhao, and X. Huang (2022) Embedding and extraction of knowledge in tree ensemble classifiers. Machine Learning 111 (5), p. 1925–1958. External Links: ISSN 1573-0565, Document Cited by: §6. C. Igel, N. Hansen, and S. Roth (2007) Covariance matrix adaptation for multi-objective optimization. Evolutionary computation 15 (1), p. 1–28. Cited by: §5.1, §5.2, §5.3. N. Jansen, F. Corzilius, M. Volk, R. Wimmer, E. Abraham, J. Katoen, and B. Becker (2014) Accelerating parametric probabilistic verification. In Quantitative Evaluation of Systems, Cham, p. 404–420. External Links: ISBN 978-3-319-10696-0 Cited by: §1, §2, §4.2. S. K. Jha, E. M. Clarke, C. J. Langmead, A. Legay, A. Platzer, and P. Zuliani (2009) A bayesian approach to model checking biological systems. In 7th Int. Conf. on Computational Methods in Systems Biology, p. 218–234. Cited by: §1, §1. T. Jin, Q. Wang, L. Xu, C. Pan, L. Dou, H. Qian, L. He, and T. Xie (2019) FinExpert: domain-specific test generation for fintech systems. In Proc. of the 27th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering, p. 853–862. Cited by: §6. D. Kerrigan, J. Hullman, and E. Bertini (2021) A survey of domain knowledge elicitation in applied machine learning. Multimodal Technologies and Interaction 5 (12), p. 73. Cited by: §6. S. Kullback and R. A. Leibler (1951) On Information and Sufficiency. The Annals of Mathematical Statistics 22 (1), p. 79 – 86. Cited by: §4.3. M. Kwiatkowska, G. Norman, and D. Parker (2004) Modelling and verification of probabilistic systems. Mathematical Techniques for Analyzing Concurrent and Probabilistic Systems. CRM Monograph Series 23, p. 93–215. Cited by: §1. M. Kwiatkowska, G. Norman, and D. Parker (2007) Stochastic model checking. Formal Methods for Performance Evaluation: 7th International School on Formal Methods for the Design of Computer, Communication, and Software Systems, SFM 2007, Bertinoro, Italy, May 28-June 2, 2007, Advanced Lectures 7, p. 220–270. Cited by: §1, §1, §2. M. Kwiatkowska, G. Norman, and D. Parker (2011) PRISM 4.0: Verification of probabilistic real-time systems. In Computer Aided Verification, G. Gopalakrishnan and S. Qadeer (Eds.), LNCS, Vol. 6806, Berlin, Heidelberg, p. 585–591. External Links: ISBN 978-3-642-22110-1 Cited by: §1, §2. M. Kwiatkowska, G. Norman, and D. Parker (2022) Probabilistic Model Checking and Autonomy. Annual Review of Control, Robotics, and Autonomous Systems 5 (1), p. 385–410. External Links: Document Cited by: §1, §6, §6. B. Lacerda, F. Faruq, D. Parker, and N. Hawes (2019) Probabilistic planning with formal performance guarantees for mobile service robots. The International Journal of Robotics Research 38 (9), p. 1098–1123. Cited by: §6. M. Lahijanian, S. B. Andersson, and C. Belta (2011) Temporal logic motion planning and control with probabilistic satisfaction guarantees. IEEE Transactions on Robotics 28 (2), p. 396–409. Cited by: §6. M. Lahijanian, S. B. Andersson, and C. Belta (2015) Formal verification and synthesis for discrete-time stochastic systems. IEEE Transactions on Automatic Control 60 (8), p. 2031–2045. Cited by: §6. N. Li, S. Adepu, E. Kang, and D. Garlan (2020) Explanations for human-on-the-loop: a probabilistic model checking approach. In Proceedings of the IEEE/ACM 15th International Symposium on Software Engineering for Adaptive and Self-Managing Systems, p. 181–187. Cited by: §1. J. C. Moreno, A. Lopes, D. Garlan, and B. Schmerl (2014) Impact models for architecture-based self-adaptive systems. In Int. Conference on Formal Aspects of Component Software, p. 89–107. Cited by: §1. V. Nardone, A. Santone, M. Tipaldi, and L. Glielmo (2016) Probabilistic model checking applied to autonomous spacecraft reconfiguration. In 2016 IEEE Metrology for Aerospace (MetroAeroSpace), p. 556–560. Cited by: §6. B. Nuseibeh and S. Easterbrook (2000) Requirements engineering: a roadmap. In Proc. of the Conference on the Future of Software Engineering, p. 35–46. Cited by: §6. A. O’Hagan, C. E. Buck, A. Daneshkhah, J. R. Eiser, P. H. Garthwaite, D. J. Jenkinson, J. E. Oakley, and T. Rakow (2006) Uncertain judgements: eliciting experts’ probabilities. John Wiley & Sons. External Links: ISBN 0-470-02999-4 Cited by: §6. G. Oghabi, J. Bentahar, and A. Benharref (2011) On the verification of behavioral and probabilistic web services using transformation. In IEEE Int. Conf. on Web Services, p. 548–555. Cited by: §6. V. Pareto (1964) Cours d’économie politique. Vol. 1, Librairie Droz. Cited by: §4.3. D. Perez-Palacin and J. Merseguer (2010) Performance evaluation of self-reconfigurable service-oriented software with stochastic petri nets. Electronic Notes in Theoretical Computer Science 261, p. 181–201. Cited by: §1. K. Sastry, D. Goldberg, and G. Kendall (2005) Genetic algorithms. In Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques, E. K. Burke and G. Kendall (Eds.), p. 97–125. Cited by: footnote 4. S. Sonnentag (1998) Expertise in professional software design: a process study.. Journal of applied psychology 83 (5), p. 703. Cited by: §6. R. Stewart and S. Ermon (2017) Label-free supervision of neural networks with physics and domain knowledge. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 31. Cited by: §6. G. Su and D. S. Rosenblum (2013) Asymptotic bounds for quantitative verification of perturbed probabilistic systems. In International Conference on Formal Engineering Methods, p. 297–312. Cited by: §6. S. Tigane, L. Kahloul, N. Hamani, M. Khalgui, and M. A. Ali (2022) On quantitative properties preservation in reconfigurable generalized stochastic petri nets. IEEE Transactions on Systems, Man, and Cybernetics: Systems. Cited by: §1. M. Tomy, B. Lacerda, N. Hawes, and J. L. Wyatt (2020) Battery charge scheduling in long-life autonomous mobile robots via multi-objective decision making under uncertainty. Robotics and Autonomous Systems 133, p. 103629. Cited by: §6. G. Vázquez, A. Evangelidis, S. Shahbeigi, R. Calinescu, and S. Gerasimou (2026) Mind the prompt: self-adaptive generation of task plan explanations via llms. In Proceedings of the 21st International Conference on Software Engineering for Adaptive and Self-Managing Systems, p. 162–175. Cited by: §6. G. Vázquez, A. Evangelidis, S. Shahbeigi, and S. Gerasimou (2025) Adaptive human-robot collaborative missions using hybrid task planning. In 2025 IEEE/ACM 20th Symposium on Software Engineering for Adaptive and Self-Managing Systems (SEAMS), p. 73–84. Cited by: §6. [82] T. Vogel Self-adaptive systems exemplars. Note: https://w.hpi.uni-potsdam.de/giese/public/selfadapt/exemplars Cited by: §6. G. Walter and T. Augustin (2009) Imprecision and prior-data conflict in generalized Bayesian inference. Journal of Statistical Theory and Practice 3 (1), p. 255–271. Cited by: §2. D. Weyns, R. Calinescu, R. Mirandola, K. Tei, M. Acosta, N. Bencomo, A. Bennaceur, N. Boltz, T. Bures, J. Camara, et al. (2023) Towards a research agenda for understanding and managing uncertainty in self-adaptive systems. ACM SIGSOFT Software Engineering Notes 48 (4), p. 20–36. Cited by: §1. G. Wright and P. Ayton (1987) Eliciting and modelling expert knowledge. Decision Support Systems 3 (1), p. 13–26. External Links: ISSN 0167-9236 Cited by: §6. X. Xie, J. Niu, X. Liu, Z. Chen, S. Tang, and S. Yu (2021) A survey on incorporating domain knowledge into deep learning for medical image analysis. Medical Image Analysis 69, p. 101985. Cited by: §1, §6. F. Zervoudakis, D. S. Rosenblum, S. Elbaum, and A. Finkelstein (2013) Cascading verification: an integrated method for domain-specific model checking. In Proc. of the 9th Joint Meeting on Foundations of Software Engineering, p. 400–410. Cited by: §6. X. Zhao, R. Calinescu, S. Gerasimou, V. Robu, and D. Flynn (2020a) Interval Change-Point Detection for Runtime Probabilistic Model Checking. In Proc. of the 35th IEEE/ACM Int. Conf. on Automated Software Engineering, ASE’20, p. 163–174. External Links: Document Cited by: §4.1, §6. X. Zhao, S. Gerasimou, R. Calinescu, C. Imrie, V. Robu, and D. Flynn (2024) Bayesian learning for the robust verification of autonomous robots. Communications Engineering 3 (1), p. 18. Cited by: §1, §5.3, §6, §6, §7. X. Zhao, M. Osborne, J. Lantair, V. Robu, D. Flynn, X. Huang, M. Fisher, F. Papacchini, and A. Ferrando (2019a) Towards integrating formal verification of autonomous robots with battery prognostics and health management. In Software Engineering and Formal Methods, P. C. Ölveczky and G. Salaün (Eds.), LNCS, Vol. 11724, p. 105–124. Cited by: §1, §6. X. Zhao, V. Robu, D. Flynn, F. Dinmohammadi, M. Fisher, and M. Webster (2019b) Probabilistic model checking of robots deployed in extreme environments. In Proc. of the 33rd AAAI Conference on Artificial Intelligence, Vol. 33, Honolulu, Hawaii, USA, p. 8076–8084. Cited by: §6, §6. X. Zhao, K. Salako, L. Strigini, V. Robu, and D. Flynn (2020b) Assessing safety-critical systems from operational testing: A study on autonomous vehicles. Information and Software Technology 128, p. 106393. External Links: ISSN 0950-5849 Cited by: §6. E. Zitzler, J. Knowles, and L. Thiele (2008) Quality assessment of Pareto set approximations. Multiobjective optimization, p. 373–404. Cited by: §5.2, §5.2. E. Zitzler, M. Laumanns, and L. Thiele (2001) SPEA2: improving the strength Pareto evolutionary algorithm. In Evolutionary Methods for Design Optimization and Control with Applications to Industrial Problems (EUROGEN’01), p. 95–100. Cited by: §5.1, §5.2, §5.3.