A MATLAB Toolbox for Fuzzy Relational Calculus in Variety of Fuzzy Algebras

1. Introduction

We present theory and software for fuzzy relational calculus (FRC) with applications in variety of areas. Attention is paid on the inverse problem resolution and obtaining complete solution set for fuzzy linear system of equations (FLSEs) and fuzzy linear system of inequalities (FLSIs). Applications are given for finite fuzzy machines – behaviour, reduction and minimization, as well as for optimization of linear objective function under FLSEs or FLSIs constraints.

Chronologically the main problem - how to solve FLSEs arose when we investigated finite fuzzy machines: in order to find finite behaviour matrix, to reduce or minimize states, we needed a solution method for FLSEs. This motivated our interest to FLSEs.

Publications for FLSEs demonstrate long period of investigations for discovering methods and procedures to solve them. Traditional linear algebra methods [1] cannot be used here – operations in fuzzy algebras are different from classical addition and multiplication.

A short review of the main results and approaches follows. Fuzzy relational equations (FREs) were first proposed in connection with medical diagnosis by Sanchez [2,3]. In [4] (1976) he gave formulas how to compute the maximal solution in the case of $max-min$ composition (Theorem 5) and the minimal solution in the case of $min-max$ composition (Theorem 6).

After Sanchez results attention was paid to the lower solutions of $max-min$ FLSEs, the complete solution set and the solution procedures. Higashi and Klir in [5] (1984) improved the results from Czogała et. all [6] (1982) and derived several schemes for solving finite FREs. assuming that the coefficients on the right-hind side of FLSEs are ordered decreasingly. The procedures are tremendous. Analogical approach was given by Miyakoshi and Shimbo [7] (1986), the results are valid for FREs with triangular norms. Pappis and Sugeno [8] (1985) proposed analytical method for obtaining the lower solutions, marking the essential coefficients for the solution and introducing the dominance condition (that helps to exclude from investigation some non-essential paths); nevertheless, the set of potential lower solutions has to be investigated in order to find valid one. Pappis and Adamopoulos [9,10,11], proposed procedures for inverse problem resolution in ZBasic for matrices with dimensions up to 10x10. Peeva [12] (1992) proposed analytical method and algorithms to solve FLSEs, introducing three semantic types of coefficients (G and E as essential and S as non-essential) for the potential solutions. She also marked that the time complexity of the algorithm for finding lower solutions is exponential, later proved by Chen and Wang [13] (2002).

After these early years the interest was paid to other compositions: $min-max$ (Peeva and Kyosev [15]); $max-$product (Bourke and Fisher [16], Di Nola et. all [17], Loetamonphong and Fang [18], Markovskii [19]); $sup-t-$norm (Bartl and Belohlávek [20], Feng Sun [21], Nosková and Perfilieva [22]); $inf-$residuum (Perfilieva and Nosková [23]), etc.

There exist many other valuable publications and generalizations for the inverse problem resolution: see Bartl as well as Bartl and co-authors [24–26]; De Baets [27]; Di Nola et. all [17]; Li and Fang [28]; Molai and Khorram [29]; Shieh [30]; Shivanian [31]; Wu and Guu [32].

Finding analytical method for the lower solutions is a long and arduous period [33], not to mention the even more difficult task of a software. The procedures are tremendous. The time complexity of the problem is exponential [13]. Various approaches were proposed to find solutions (extremal and the set of all solutions):

• analytical methods based on algebraic-logical approach [5,7,8,14,15,17,24,27,28], with applications in various subjects [15] and with the first software [34] by Kyosev;
• characteristic matrix and decomposition [7,8];
• covering and binding variables [19,35–37];
• partitions and irreducible paths [5];
• solution based matrix [13], etc.

In this paper attention is paid on software description (based on our unified approach and exact methods) for solving FLSEs in some$BL-$algebras: in Gödel algebra, in Goguen algebra and in Łukasiewicz algebra.

In Section 2 we give the basic notions for the exposition - $BL-$algebra, Gödel algebra, Goguen algebra and Łukasiewicz algebra, $t-$ and $s-$ norms, compositions of fuzzy relations and variety of products for finite fuzzy matrices.

Section 3 presents the theory and exact method for solving FLSEs or FLSIs in some $BL-$algebras. A short overview is given of the chronological development of the fundamental ideas.

Applications for finite fuzzy machines and for linear optimization problems are subject of Section 4. They concern solving behaviour, reduction and minimization problems for finite fuzzy machines and solving linear optimization problems subject to FLSEs of FLSIs constraints.

The exact algorithms currently used in the software [38] are explained in details in Section 5, together with their computational complexity.

Section 6 presents the software from [38] itself and describes each of its modules and capabilities.

All algorithms and software functions are supported by examples provided in Appendix A.

2. Basic Notions

In this section the notions for classical algebra, orders and lattices are given according to [1,39], for $BL-$algebra – according to [40], for fuzziness - according to [15].

Partial order relation on a partially ordered set (poset) P is denoted by the symbol ≤. By a greatest element of a poset P we mean an element $b\in P$ such that $x\le b$ for all $x\in P$. The least element of P is defined dually.

2.1. $BL-$algebra

$BL-$algebra is the algebraic structure:

$$BL=〈L,\vee ,\wedge ,\ast ,\to ,0,1〉,$$

where $\vee ,\wedge ,\ast ,\to$ are binary operations, $0,1$ are constants and:

• $L=〈L,\vee ,\wedge ,0,1〉$ is a lattice with universal bounds 0 and 1;
• $L=〈L,\ast ,1〉$ is a commutative semigroup;
• * and → establish an adjoint couple:

  $$z\le \left(x\to y\right)\iff x\ast z\le y,\forall x,y,z\in L.$$
• for all $x,y\in L$

  $$x\ast (x\to y)=x\wedge yand(x\to y)\vee (y\to x)=1.$$

We suppose in next exposition that $L=[0,1]$ and $x,y\in [0,1]$.

The following algebras are examples for $BL-$algebras:

• Gödel algebra

  $$B{L}_G=〈[0,1],\vee ,\wedge ,{\to }_G,0,1〉,$$

  where operations are

  (a)

  Maximum or $s_3-$conorm:

  $$max\{x,y\}=x\vee y.$$

  (1)

  (b)

  Minimum or $t_3-$norm:

  $$min\{x,y\}=x\wedge y.$$

  (2)

  (c)

  The residuum ${\to }_G$ is

  $$x{\to }_{G}y=\left\{\begin{array}{c}1ifx\le y\hfill \\ yifx>y\hfill \end{array}\right..$$

  (3)

  (d)

  A supplementary operation is useful

  $$x\epsilon y=\left\{\begin{array}{c}y,ifx<y\hfill \\ 0,ifx\ge y\hfill \end{array}\right..$$

  (4)

• Product (Goguen) algebra

  $$B{L}_P=〈[0,1],\vee ,\wedge ,\circ ,{\to }_P,0,1〉,$$

  where max and min are as (1) and (2), respectively, ∘ is the conventional real number multiplication (the $t_2-$ norm, i. e, $t_2(x,y)=xy$) and the residuum ${\to }_P$ is

  $$x{\to }_{P}y=\left\{\begin{array}{c}1ifx\le y\hfill \\ {\displaystyle \frac{y}{x}}ifx>y\hfill \end{array}\right..$$

  (5)
  Here the supplementary useful operation is:

  $$x\gamma y=\left\{\begin{array}{c}0ifx\ge y\hfill \\ {\displaystyle \frac{y-x}{1-x}}ifx<y\hfill \end{array}\right..$$

  (6)
• ukasiewicz algebra

  $$B{L}_L=〈[0,1],\vee ,\wedge ,\otimes ,{\to }_L,0,1〉,$$

  where max and min are as (1) and (2), respectively, and

  (a)

  $x\otimes y=0\vee (x+y-1)\equiv t_1(x,y).$

  (b)

  The residuum ${\to }_L$ is

  $$x{\to }_{L}y=1\wedge (1-x+y).$$

  (7)

  (c)

  A supplementary operation is useful

  $$x\delta y=0\vee (y-x).$$

  (8)

These operations, as well as the operations in Table 1 are realized in [38], see Section 5, Section 6, and Apendix Appendix A for details and working examples.

These three couples are also described as $t-$norms and $s-$norms (or $t-$conorms):

2.2. Compositions of Fuzzy Relations and Fuzzy Matrix Products

Let X and Y be crisp sets and $F(X\times Y)$ denote all fuzzy sets over $X\times Y$. A fuzzy relation$R\in F(X\times Y)$ is defined as a fuzzy subset of the Cartesian product $X\times Y$,

$$R=\left\{\left((x,y),{\mu }_R(x,y)\right)(x,y)\in X\times Y,{\mu }_R:X\times Y\to [0,1]\right\}.$$

The inverse(or transpose) $R^{-1}\in F(Y\times X)$ of $R\in F(X\times Y)$ is defined as

$$R^{-1}(y,x)=R(x,y)forallpairs(y,x)\in Y\times X.$$

The finite fuzzy matrix (FFM) $A={\left(a_{ij}\right)}_{m\times n}$, with $a_{ij}\in [0,1]$ for each $i,j,1\le i\le m$, $1\le j\le n$, is called a membership matrix.

Any fuzzy relation $R\in F(X\times Y)$ over finite support $X\times Y$ is representable by FFM, written for convenience with the same letter $R={\left(r_{ij}\right)}_{m\times n}$, where $r_{ij}={\mu }_R(x_i,y_j)$ for any $(x_i,y_j)\in X\times Y$.

The matrix representation of the inverse fuzzy relation is $R^{-1}=\left(r_{ji}\right)$, called transpose or inverse of the given matrix. In this sense $R^{-1}=R^t$.

Two FFMs are called conformable, if the number of the columns in the first FFM equals the number of the rows in the second FFM: the matrices $A={\left(a_{ij}\right)}_{m\times p}$ and $B={\left(b_{ij}\right)}_{p\times n}$ are conformable and their product $C={\left(c_{ij}\right)}_{m\times n}=AB$, in this order, makes sense.

We consider some products of FFMs in $BL-$algebras. Since in references special attention is paid on $max-min$, $min-max$ and $max-product$ compositions, they are specially listed in next Definition 1.

Definition 1.

Let $A={\left(a_{ij}\right)}_{m\times p}$ and $B={\left(b_{ij}\right)}_{p\times n}$ be conformable FFMs.

The FFM $C={\left(c_{ij}\right)}_{m\times n}$, where

i) 

$C=A•B,$ is called the $max-min$ product of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\vee }}(a_{ik}\wedge b_{kj})when1\le i\le m,1\le j\le n.$$

ii) 

$C=A\circ B,$ is called the $min-max$ product of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\wedge }}(a_{ik}\vee b_{kj})when1\le i\le m,1\le j\le n.$$

iii) 

$C=A\odot B,$ is called the $max-$product (the product ’.’ is the conventional real number multiplication) of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\vee }}(a_{ik}.b_{kj})when1\le i\le m,1\le j\le n.$$

iv) 

$C=A{\to }_{G}B,$ is called the $min-{\to }_G$ product of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\wedge }}\left(a_{ik}{\to }_G{b}_{kj}\right)when1\le i\le m,1\le j\le n.$$

v) 

$C=A{\to }_{P}B,$ is called the $min-{\to }_P$ product of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\wedge }}\left(a_{ik}{\to }_P{b}_{kj}\right)when1\le i\le m,1\le j\le n.$$

vi) 

$C=A\epsilon B,$ is called the $max-\epsilon$ product of A and B if

$$c_{ij}=\underset{k=1}{\overset{p}{\vee }}\left(a_{ik}\epsilon b_{kj}\right)when1\le i\le m,1\le j\le n.$$

For simplicity and compactness of exposition we write $AB$ for any of the matrix products introduced in Definition 1.

Definition 2.

Let $A={\left({\mu }_{ij}^A\right)}_{m\times p}$ and $B={\left({\mu }_{ij}^B\right)}_{p\times n}$ be conformable FFMs. The matrix $C={\left({\mu }_{ij}^C\right)}_{m\times n}$ is called$s-t$ product of A and B in $BL-$algebra, written $C=A{\ast }_{BL}B$, if for each $i,j$, $1\le i\le m,1\le j\le n$ it holds:

$${\mu }_{ij}^C={s_r}_{k=1}^p\left(t_r({\mu }_{ik}^A,{\mu }_{kj}^B)\right),$$

(9)

where $s_r$ and $t_r$ are $s-$norm and $t-$norm respectively for $r=1,2,3$, see Table 1.

Definition 3.

Let $A={\left({\mu }_{ij}^A\right)}_{m\times p}$ and $B={\left({\mu }_{ij}^B\right)}_{p\times n}$ be conformable FFMs. The matrix $C={\left({\mu }_{ij}^C\right)}_{m\times n}$ is called $t-{\to }_{BL}$ product of A and B in $BL-$algebra, written $C=A{\to }_{BL}B$, if for each $i,j$, $1\le i\le m,1\le j\le n$ it holds:

$${\mu }_{ij}^C={t_3}_{k=1}^p\left(\left({\mu }_{ik}^A{\stackrel{r}{\to }}_{BL}{\mu }_{kj}^B\right)\right),$$

(10)

where $t_3$ is the $t-$norm (see Table 1) and ${\stackrel{r}{\to }}_{BL}=\delta$ for $r=1$, see (8); ${\stackrel{r}{\to }}_{BL}=\gamma$ for $r=2$, see (6); ${\stackrel{r}{\to }}_{BL}=\epsilon$ for $r=3,$ see (4).

2.3. Direct and inverse problem resolution

Let $A={\left(a_{ij}\right)}_{m\times p}$ and $B={\left(b_{ij}\right)}_{p\times n}$ be conformable FFMs. Computing their product according to Definition 1, 2 or 3, respectively, is called direct problem resolution.

Let Let $A={\left(a_{ij}\right)}_{m\times p}$ and $C={\left(c_{ij}\right)}_{m\times n}$ be given FFMs. Computing the matrix $B={\left(b_{ij}\right)}_{p\times n}$, such that $AB=C$ or $A{\ast }_{BL}B=C$ or $A{\to }_{BL}B=C$ is called inverse problem resolution.

Details for algorithm and software on this subject are given in subSection 5.1 and Section 6.2 with working examples in Appendix A.1.

3. FLSEs and FLSIs in Some $BL-$Algebras

We began our investigations before 1992 [12] with $max-min$, i.e. $s_3-t_3$-norm FLSEs, given with long form and short form description:

$$\left|\begin{array}{ccccc}(a_{11}\wedge x_1)& \vee \dots \vee & (a_{1n}\wedge x_n)& =& b_1\\ \dots & \dots & \dots & \dots & \dots \\ (a_{m1}\wedge x_1)& \vee \dots \vee & (a_{mn}\wedge x_n)& =& b_m\end{array}\leftrightarrow A•X=B\right..$$

(11)

Here $A={\left(a_{ij}\right)}_{m\times n}$ stands for the matrix of coefficients, $X={\left(x_j\right)}_{n\times 1}$ – for the matrix of unknowns, $B={\left(b_i\right)}_{m\times 1}$ is the right-hand side of the system, $a_{ij}$, $b_i$, $x_j\in [0,1]$ for $i=1,...,m$ and $j=1,...,n.$

The next step was to investigate $min-max$ FLSEs $A\circ X=B$ and $max-product$ FLSEs $A\odot X=B$. The main author’s results for all these FLSEs are published in [41].

The theoretical background - how to solve (11), is implemented as universal scheme in theory, algorithms and software for all the other kind of compositions for the FLSEs and FLSIs ($s_3-t_3$, $s_3-t_2$, $s_3-t_1$ as well as $t_3-s_3$, $t_3-s_2$, $t_3-s_1$), presented in the next subsection.

3.1. $s_3-t$-norm FLSEs in $BL-$algebras

The natural extension of the theory for $max-min$ FLSEs (11) leads to the unified approach for FLSEs when the composition is $s_3-t$:

$$\left|\begin{array}{c}{t}_r(a_{11},x_1)\vee t_r(a_{12},x_2)\vee \dots \vee t_r(a_{1n},x_n)=b_1\hfill \\ t_r(a_{21},x_1)\vee t_r(a_{22},x_2)\vee \dots \vee t_r(a_{2n},x_n)=b_2\hfill \\ \dots \hfill \\ t_r(a_{m1},x_1)\vee t_r(a_{m2},x_2)\vee \dots \vee t_r(a_{mn},x_n)=b_m\hfill \end{array}\leftrightarrow A{\ast }_{BL}X=B\right..$$

(12)

where $A={\left(a_{ij}\right)}_{m\times n}$ is the matrix of coefficients, $B={\left(b_i\right)}_{m\times 1}$ holds for the right-hand side vector, $X={\left(x_j\right)}_{1\times n}$ is the vector of unknowns; for the indices we suppose $i=1,...,m$, $j=1,...,n$; $a_{ij},b_i\in [0,1]$ are given, $x_j\in [0,1]$ marks the unknowns in the system and $t_r$ ($r=1,2,3$) is according to Table 1.

The approach is also valid for FLSEs $A{\to }_{BL}X=B$.

Solving (12) for the unknown matrix X is called inverse problem resolution for FLSEs $A{\ast }_{BL}X=B$.

In the next exposition we suppose that the matrices A, X and B are as follows: $A={\left(a_{ij}\right)}_{m\times n}$, $X={\left(x_j\right)}_{1\times n}$ and $B={\left(b_i\right)}_{m\times 1}$.

3.2. Solutions

For FFMs $X_{p\times n}$ and $Y_{p\times n}$ the inequality

$$X\le Y$$

means ${\mu }_{ij}\left(x\right)\le {\mu }_{ij}\left(y\right)$ for each $i=1,\dots ,p$, $j=1,\dots ,n$.

Definition 4.

For the FLSEs (12) $A{\ast }_{BL}X=B$:

i) 

$X_{n\times 1}^0$ is called a point solution of $A{\ast }_{BL}X=B$ if $A{\ast }_{BL}{X}^0=B$ holds.

ii) 

The set of all solutions of (12) is called complete solution set and is denoted by $\mathbb{X}$. If $\mathbb{X}\ne \varnothing$ then (12) is called consistent, otherwise it is called inconsistent .

iii) 

A solution $\widehat{X}\in \mathbb{X}$ is called an upper or maximal solution of (12) if for any $X\in \mathbb{X}$ the relation $\widehat{X}\le X$ implies $X=\widehat{X}$. When the upper solution is unique, it is called the greatest or maximumsolution .

iv) 

A solution $\check{X}\in \mathbb{X}$ is called a lower or minimal solution of (12) if for any $X\in \mathbb{X}$ the relation $X\le \check{X}$ implies $X=\check{X}$. When the lower solution is unique, it is called the least or minimum solution .

v) 

$(X_1,\dots ,X_n)$ with $X_j\subseteq [0,1]$ for each j, $1\le j\le n$, is called an interval solution if any $X^0={\left(x_j^0\right)}_{n\times 1}$ belongs to $\mathbb{X}$ when $x_j^0\in X_j$ for each j, $1\le j\le n$.

vi) 

Any interval solution, whose interval bounds are bounded by a lower solution from the left and by the greatest solution from the right, is called maximal interval solution .

3.3. Solvability and greatest solution

We give the analytical expression for the greatest solution of FLSEs (12) and a criterion for its consistency.

Theorem 1.

[14] Let A and B be FFMs and let $\mathbb{X}$ be the complete solution set of $A{\ast }_{BL}X=B$. Then:

i) 

$\mathbb{X}\ne \varnothing \iff {\widehat{X}}_{BL}=A^t{\to }_{BL}B\in \mathbb{X}$;

ii) 

If the FLSEs (12) is consistent, then ${\widehat{X}}_{BL}=A^t{\to }_{BL}B$ is its greatest solution.

iii) 

There exists polynomial time algorithm for computing $A^t{\to }_{BL}B$.

Corollary 1.

[14] The following statements are valid for FLSEs $A{\ast }_{BL}X=B$:

i) 

The FLSEs (12) is consistent iff $B=A{\ast }_{BL}\left(A^t{\to }_{BL}B\right)$ holds;

ii) 

There exists polynomial time algorithm for establishing consistency of the FLSEs (12) and for computing its greatest solution ${\widehat{X}}_{BL}=A^t{\to }_{BL}B.$

In fact, Theorem 1, item ii), gives the analytical expression for the greatest solution of the equation (12). It is the core of the algorithms and software presented in [38] for establishing consistency of FLSEs and finding its extremal solution: maximum for the equation (12) and minimum for the dual case.

First we present finding the greatest solution of FLSEs (12) $A{\ast }_{BL}X=B$ , a criterion for its consistency and and how to find the lower solutions, based on [42].

3.4. Greatest Solution

According to Theorem 1 any solvable $max-t$-norm FLSEs has unique greatest solution. In order to find all solutions of the consistent system, it is necessary to find its greatest solution and all of its minimal solutions. Finding the greatest solution is relatively simple task which can be used as a criteria for establishing consistency of the system. Finding all minimal solutions is reasonable for consistent system.

3.4.1. Classical Approach

The approach how to solve (12) is based on Theorem 1 and Corollary 1.

If (12) is consistent, its greatest solution is $\widehat{X}=\left({\widehat{x}}_j\right)=A^t{\to }_{BL}B$.

Using this fact, an appropriate algorithm for checking consistency of the system and for finding its greatest solution can be obtained. Its computational complexity is $O(m.n^2)$. Nevertheless that it is simple, it is too hard for such a task.

3.4.2. More Efficient Approach

Here we propose a simpler way to answer both questions for (12) – consistency and the greatest solution (simultaneously computing the greatest solution and establishing consistency of (12)). Extending the ideas of [12] and [15], we work with four types of coefficients (S, E, G and H) and a boolean vector (IND).

In the FLSEs (12):

• If $r=3$ i.e. operation is minimum ($t_3(x,y)=min\{x,y\}$):

  −

  $a_{ij}$ is called S-type coefficient if $a_{ij}<b_i$.

  −

  $a_{ij}$ is called E-type coefficient if $a_{ij}=b_i$.

  −

  $a_{ij}$ is called G-type coefficient if $a_{ij}>b_i$.

  −

  $a_{ij}$ is called H-type coefficient if $a_{ij}\ge b_i$.

• If $r=2$ i.e. operation is the algebraic product ($t_2(x,y)=xy$):

  −

  $a_{ij}$ is called S-type coefficient if $a_{ij}<b_i$.

  −

  $a_{ij}$ is called E-type coefficient if $a_{ij}=b_i=0$.

  −

  $a_{ij}$ is called G-type coefficient if $a_{ij}\ge b_i>0$.

  −

  $a_{ij}$ is called H-type coefficient if $a_{ij}\ge b_i$.

• If $r=1$ i.e. operation is the ukasiewicz $t-$norm ($t_1(x,y)=max\{x+y-1,0\}$):

  −

  $a_{ij}$ is called S-type coefficient if $a_{ij}-1>b_i$.

  −

  $a_{ij}$ is called E-type coefficient if $a_{ij}=b_i=0$.

  −

  $a_{ij}$ is called G-type coefficient if $a_{ij}\ge b_i>0$.

  −

  $a_{ij}$ is called H-type coefficient if $a_{ij}-1\le b_i$.

The next Algorithm 1 uses the fact that it is possible to find the value of the unknown ${\widehat{x}}_j$ only by the $j^{th}$ column of the matrix A. For every $i=1,...,m$ and depending on the operation, we consider the following cases:

• When the operation is $t_3$

  −

  If $a_{ij}$ is E-type coefficient then the i-th equation can be satisfied by $a_{ij}\wedge x_j$ when $x_j\ge b_i$ because $a_{ij}\wedge x_j=b_i\wedge x_j=b_i$.

  −

  If $a_{ij}$ is G-type coefficient then the i-th equation can be satisfied by $a_{ij}\wedge x_j$ only when $x_j=b_i$ because $a_{ij}\wedge x_j=a_{ij}\wedge b_i=b_i$.

  −

  If $a_{ij}$ is S-type coefficient then the i-th equation cannot be satisfied by $a_{ij}\wedge x_j$ for any $x_j\in [0,1]$.

• When the operation is $t_2$

  −

  If $a_{ij}$ is E-type coefficient then the i-th equation can be satisfied by $a_{ij}{x}_j$ when $x_j\in [0,1]$ because $a_{ij}{x}_j=0x_j=b_i=0$.

  −

  If $a_{ij}$ is G-type coefficient then the i-th equation can be satisfied by $a_{ij}{x}_j$ only when $x_j=b_i/a_{ij}$ because $a_{ij}{x}_j=a_{ij}(b_i/a_{ij})=b_i$.

  −

  If $a_{ij}$ is S-type coefficient then the i-th equation cannot be satisfied by $a_{ij}{x}_j$ for any $x_j\in [0,1]$.

• When the operation is $t_1$

  −

  If $a_{ij}$ is E-type coefficient then the i-th equation can be satisfied by $a_{ij}\otimes x_j$ when $x_j\in [0,1]$ because $a_{ij}\otimes x_j=max\{0+x_j-1,0\}=max\{x_j-1,0\}=b_i=0$.

  −

  If $a_{ij}$ is H-type coefficient then the i-th equation can be satisfied by $a_{ij}\otimes x_j$ only when $x_j=1-a_{ij}+b_i$ because $a_{ij}\otimes x_j=a_{ij}+1-a_{ij}+b_i-1=b_i$.

  −

  If $a_{ij}$ is S-type coefficient then the i-th equation cannot be satisfied by $a_{ij}\otimes x_j$ for any $x_j\in [0,1]$.

Hence, S-type coefficients are not interesting because they do not lead to solution.

For the purposes of the next theorem, ${\widehat{b}}_j$ is introduced as follows:

• If the operation is $t_3$:

  $${\widehat{b}}_j=\left\{\begin{array}{c}\underset{i=1}{\overset{m}{min}}\left\{b_i\right\},\mathrm{for}\mathrm{all}\mathit{i}\mathrm{such}\mathrm{that}{a}_{ij}>b_i\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\right.$$

  (13)
• If the operation is $t_2$:

  $${\widehat{b}}_j=\left\{\begin{array}{c}\underset{i=1}{\overset{m}{min}}\{b_i/a_{ij}\},\mathrm{for}\mathrm{all}\mathit{i}\mathrm{such}\mathrm{that}{a}_{ij}>b_i\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\right.$$

  (14)
• If the operation is $t_1$:

  $${\widehat{b}}_j=\left\{\begin{array}{c}\underset{i=1}{\overset{m}{min}}\{1-a_{ij}+b_i\},\mathrm{for}\mathrm{all}\mathit{i}\mathrm{such}\mathrm{that}{a}_{ij}-1<b_i\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\right.$$

  (15)

Theorem 2.

[42] The FLSEs (12) is consistent iff $\widehat{X}=\left({\widehat{b}}_j\right)$ is its solution.

Corollary 2.

[42] In a consistent FLSEs (12):

• choosing $x_j>{\widehat{b}}_j$ for at least one $j=1,...,n$ makes the system inconsistent;
• for every $j=1,...,n$, the greatest admissible value for $x_j$ is ${\widehat{b}}_j$;
• its greatest solution is $\widehat{X}=\left({\widehat{x}}_j\right)=\left({\widehat{b}}_j\right)$, $j=1,...,n$.

In general, Theorem 2 and its Corollary show that instead of calculating $\widehat{X}=A^t{\to }_{BL}B$ we can implement a faster way to obtain $\widehat{X}=\left({\widehat{x}}_j\right)=\left({\widehat{b}}_j\right)$ (as realized in Algorithm 1).

$\widehat{X}$ is only the eventual greatest solution of the FLSEs (12) and it can be obtained for any FLSEs, even if the system is inconsistent, so the eventual solution should be checked in order to confirm that it is a solution of (12). Explicit checking for the eventual solution will increase the computational complexity of the algorithm. To avoid this in the next presented Algorithm 1 this is done by extraction of the coefficients of the potential greatest solution. For every $\left({\widehat{x}}_j\right)\in \widehat{X}$ we check which equations of (12) are satisfied (hold in the boolean vector $IND$). If at the end of the algorithm all the equations of FLSEs (12) are satisfied (i.e. all the coefficients in $IND$ are set to $TRUE$) this means that the system is consistent and the computed solution is its greatest solution, otherwise the FLSEs (12) is inconsistent.

Algorithm 1 Greatest solution and consistency of FLSEs (12).

Step 1: Initialize the vector $\widehat{X}=\left({\widehat{x}}_j\right)$ with ${\widehat{x}}_j=1$ for $j=1,...,n$. Step 2: Initialize a boolean vector $IND$ with $IN{D}_i=FALSE$ for $i=1,...,m$. This vector is used to mark equations that are satisfied by the eventual greatest solution. Step 3: For each column $j=1,...,n$ in A: walk successively through all coefficients $a_{ij},i=1,...,m$ seeking the smallest G-type coefficient. Step 4:   If $a_{ij}$ is E-type coefficient, according to (13), (14) or (15) respectively (depending on the operation), it means that the $i^{th}$ equation in the system can be solved through this coefficient, but ${\widehat{b}}_j$ still should be found. Correct $IN{D}_i$ to $TRUE$. Step 5:   For the smallest H-type (which is not E-type) coefficient correct $IN{D}_i$ to $TRUE$. All other H-type (which are not E-type) coefficients are now insignificant, as they are not leading to a solution. In $\widehat{X}$ correct the value for ${\widehat{x}}_j={\widehat{b}}_j$. Step 6: Check if all components of $IND$ are set to $TRUE$. Step 7:   If $IN{D}_i=FALSE$ for some i the system $A{\ast }_{BL}X=B$ is inconsistent. Step 8:   If $IN{D}_i=TRUE$ for all $i=1,...,m$ the system $A{\ast }_{BL}X=B$ is consistent and its greatest solution is stored in $\widehat{X}$.

Theorem 2 and its Corollary provide that if (12) is consistent, $\widehat{X}$ computed by Algorithm 1 is its greatest solution.

Algorithm 1 uses $IND$ vector to check which equations are satisfied by the eventual $\widehat{X}$. At the end of the algorithm if all components in $IND$ are $TRUE$ then $\widehat{X}$ is the greatest solution of the system, otherwise the system is inconsistent.

3.5. Lower Solutions

Each equation in (12) can be satisfied only by the terms with H-type coefficients. Also, the minimal value for every component in the solution is either the value of the corresponding $\widehat{b}$ or 0. Along these lines, the hearth of the presented here algorithm is to find H-type components $a_{ij}$ in A and to give to $X_{lo{w}_j}$ either the value of the corresponding $\widehat{b}$ when the coefficient contributes to solve the system or 0 when it doesn’t.

Using this, the set of candidates for solutions can be obtained. All candidate solutions are of three different types:

• Lower solution;
• Non-lower solution;
• Not solution at all.

The aim of the next Algorithm 2 is to extract all lower solutions and to skip the second and third types (non-lower solution or not solution at all). In order to extract all lower solutions a method, based on the idea of the dominance matrix in combination with list manipulation techniques is developed here.

3.5.1. Domination

For the purposes of presented here Algorithm 2, the definition for domination is given.

Definition 5.

Let $a_l$ and $a_k$ be the $l^{th}$ and the $k^{th}$ equations, respectively, in (12) and $b_l\ge b_k$. Equation $a_l$ is called dominant to $a_k$ and equation $a_k$ is called dominated by $a_l$, if for each $j=1,...,n$ it holds: if $a_{lj}$ is H-type coefficient then $a_{kj}$ is also H-type coefficient.

3.5.2. Extracting Lower Solutions

Lower solutions are extracted by removing from A the dominated rows. A new matrix is produced and marked with $\tilde{A}=\left(a_{\tilde{i}j}\right)$ where $\tilde{i}=1,...,\tilde{m}$, $\tilde{m}<m$ for obvious reasons. It preserves all the needed information from A to obtain the solutions.

Extraction introduced here is based on the following recursive principle. If in the $j^{th}$ column of $\tilde{A}$ there are one or more rows (${\tilde{i}}^{\ast }$) such that coefficients $a_{{\tilde{i}}^{\ast }j}$ are H-type then $x_j$ should be taken equal to the smallest $b_{{\tilde{i}}^{\ast }}$ and all rows ${\tilde{i}}^{\ast }$ should be removed from $\tilde{A}$. The same procedure is repeated for ${(j+1)}^{th}$ column of the reduced $\tilde{A}$. "Backtracking" based algorithm using this principle is presented next:

Algorithm 2 Extract the lower solutions from $\tilde{A}$.

Step 1: Initialize solution vector $X_{lo{w}_0}\left(j\right)=0,j=1,...,n$. Step 2: Initialize a vector $rows\left(\tilde{i}\right),i=1,...,\tilde{m}$ which holds all consecutive row numbers in $\tilde{A}$. This vector is used as a stopping condition for the recursion. Initially it holds all the rows in $\tilde{A}$. On every step some of the rows there are removed. When $rows$ is empty the algorithm exits from the current recursive branch. Step 3: Initialize $sols$ to be the empty set of vectors, which is supposed to be the set of all minimal solutions for current problem. Step 4: Check if $rows=\varnothing$. If so, add $X_{lo{w}_0}$ to $sols$ and exit current step. Step 5: Fix $\tilde{i}$ equal to the first element in $rows$, then for every $j=1,...,n$ such that $a_{\tilde{i}j}$ is H-type coefficient Step 6:   Create a copy of $X_{lo{w}_0}$ and update its $j^{th}$ coefficient to be equal to $b_{\tilde{i}}$. Create a copy of $rows$. Step 7:   For all $\tilde{k}$ in $rows$ if $a_{\tilde{k}j}$ is a H-type coefficient, remove $\tilde{k}$ from the copy of $rows$. Step 8:   Go to step 4 with copied in this step $rows$ and $X_{lo{w}_0}$, i.e. start new recursive branch with reduced $rows$ and changed $X_{lo{w}_0}$.

3.6. General algorithm

The next algorithm is based on the above given Algorithms 1 and 2.

Algorithm 3 Solving $A{\ast }_{BL}X=B$.

Step 1: Obtain input data for the matrices A and B. Step 2: Obtain the greatest solution for the system and check it for consistency (Algorithm 1). Step 3: If the system is inconsistent - Exit. Step 4: Obtain the matrix $\tilde{A}$. Step 5: Obtain all minimal solutions from $\tilde{A}$ and B, Algorithm 2.

Algorithm 2 (Step 5) is the slowest part of the Algorithm 3. In general this algorithm has its best and worst cases and this is the most important improvement according to algebraic-logical approach [41] (from the time complexity point of view). Algorithm 2 is going to have the same time complexity as the algorithm presented in [14] only in its worst case.

Concerning the FLSEs $A{\to }_{BL}X=B$ - this is the analogical case and we do not describe the theory here, but software is also available in [38].

More details for algorithms and software for all of the main problems in this section are given in subSection 5.2 and Section 6.3 with working examples in Apendix Appendix A.2.

Software for inverse problem resolution is provided in [34] (only for $max-min$ and $min-max$ FLSEs) and in [38] for all other cases.

3.7. $s_3-t$-Norm FLSIs in $BL-$Algebras

For the FLSIs $A{\ast }_{BL}X\le B$, $A{\ast }_{BL}X\ge B$, $A{\to }_{BL}X\le B$, $A{\to }_{BL}X\ge B$, we implement analogical approach as described for the FLSEs $A{\ast }_{BL}=B$.

For algorithms and software for FLSIs see also subSection 5.4, Section 5.5 and Section 6.3 with working examples in Apendix Appendix A.2.

The algorithms implemented in the software [38] are described in Section 5. They are the result of more than a decade of continuous refinement [14,15]. Earlier versions were built on the analytical methods [12,14,15], aimed at solving FLSEs, especially for the max-min composition.

3.8. Short Chronological Overview

The historical core of these methods was the algebraic-logical approach [12,14,15,41] which is still the main basis of all the currently used algorithms. This approach introduced a structured methodology involving potential extremal solutions, consistency checks via index vectors, and coefficient classification. The coefficients in the matrix of the coefficients A were categorized into three semantic types : S-type, E-type, and G-type coefficients [12,15]. A later refinement introduced H-type coefficients [41], which acted as a lightweight abstraction to further streamline certain matrix operations without changing the semantic interpretation of the system.

Such classifications made it possible to construct associated matrix [12] that was used to generate solution intervals and identify dominating paths. The dominance relation [41] was fundamental for reducing complexity by pruning branches that could not lead to solutions.

The next fundamental step is the construction of so-called help matrix [15], used to guide the recursive generation of candidate solutions. In some algorithmic versions, these structures were complemented by list-based procedures, where possible solutions and operations were represented and manipulated as lists rather than fixed-size matrices. These list-based algorithms enabled greater flexibility in traversing the solutions, but turn out to be even more complex from implementational point of view. Although the computational complexity of this problem is exponential [13], various simplification heuristics were developed to make the process more tractable in practice.

Over time, these methods evolved and expanded to support a wide variety of compositions in BL-algebras, including not only max–min but also $min-max$, max-product, and other BL systems such as those using the Goguen or Łukasiewicz implications. [42–44]. New versions of the algorithms integrated modular procedures that preserved the core logical structure while allowing extensibility to different algebras.

It is worth mentioning that the current implementation still employs a structure referred to as the help matrix, although its construction and role have significantly changed. The current algorithmic flow is explained in Section 5.

The evolution of these methods and their application to broader classes of fuzzy systems is documented in a series of publications [42–44]. These works include concrete examples and software implementations. Although the internal architecture of the software has shifted toward more unified and abstract representations (as detailed in Section 6), many of the core ideas — such as the use of coefficient classification, list-based construction, and filtering - continue to influence the current implementation.

4. Some Applications

4.1. Linear Dependency/Independency

Let $A\left(1\right)={\left(a_{i1}\right)}_{m\times 1}$, $A\left(2\right)={\left(a_{i2}\right)}_{m\times 1},\dots ,A\left(n\right)={\left(a_{in}\right)}_{m\times 1}$ be fuzzy column-vectors with elements $a_{ij}\in [0,1]$ for $i=1,...,m$, $j=1,...,n$.

The fuzzy column-vector $B={\left(b_j\right)}_{m\times 1}$ is called linear combination of $A\left(1\right),A\left(2\right),\dots ,A\left(n\right)$ if there exist $X={\left(x_j\right)}_{1\times n}$ such that $A{\ast }_{BL}X=B.$ Checking whether $B_{m\times 1}$ is a linear combination of $A\left(1\right),$$\dots ,A\left(n\right)$ requires to solve the system (12)

$$A{\ast }_{BL}X=B$$

for the unknown X, if $A_{m\times n}$ has as columns $A\left(1\right)$, $A\left(2\right),$$\dots ,$$A\left(n\right).$

If the FLSEs $A{\ast }_{BL}X=B$ is consistent the vectors $A\left(1\right),A\left(2\right),\dots ,A\left(n\right),B$ are linearly dependent. When the FLSEs $A{\ast }_{BL}X=B$ is inconsistent, the vectors $A\left(1\right),A\left(2\right),\dots ,A\left(n\right),B$ are linearly independent.

We implement this for instance in subSection 4.2 for finite fuzzy machines when computing the finite behaviour matrix, as well as for minimization of states.

The algorithms in subSection 5.6 and Section 5.7 are for establishing linear combination as well as linear dependency/independency. Software is described in subSection 6.2. Examples are given in Appendix A.1.

4.2. Finite Fuzzy Machines

We present finite fuzzy machines (FFMachs) with emphasis on the complete behavior matrix, reduction and minimization problems. For solving all these problems we implement the software from [38],

FFMachs were introduced by Santos, see [45] – 47]. Here the presentation is based on our results as published in [15,48].

For a finite set C we denote by $\left|C\right|$ its cardinality.

Definition 6.

A finite fuzzy machine is a quadruple

$$\mathcal{A}=(X,Q,Y,\mathcal{M}),$$

where:

i) 

$X,Q,Y$ are nonempty finite sets of input letters, states and output letters, respectively.

ii) 

$\mathcal{M}$ is the set of transition-output matrices of $\mathcal{A}$, that determines its stepwise behavior. Each matrix $M\left(x\right|y)=\left(m_{q{q}^{\prime }}\left(x\right|y)\right)\in \mathcal{M}$ is a square matrix of order $\left|Q\right|$ and $x\in X,y\in Y,q,q^{\prime }\in Q,m_{q{q}^{\prime }}\left(x\right|y)\in [0,1].$

In Definition 6 ii) we mark by $\left(x\right|y)$ the pair $(x,y)\in X\times Y$ to emphasize that $x\in X$ is the input letter when $y\in Y$ is the output response letter. The same stipulation is used in next exposition for input-output pair of words of the same length.

We regard $m_{q{q}^{\prime }}\left(x\right|y)$ as the degree of membership for the FFMach to enter state $q^{\prime }\in Q$ and produce output $y\in Y$ if the present state is $q\in Q$ and the input is $x\in X$.

4.2.1. Extended input-output behavior

While $\mathcal{M}$ is the set of transition-output matrices that describes operating of $\mathcal{A}$ for exactly one step, in this subsection we will be interested in matrices that describe operating of $\mathcal{A}$ for input-output pairs of words.

The free monoid of the words over the finite set X is denoted by $X^{\ast }$ with the empty worde as the identity element. If $X\ne \varnothing$ then $X^{\ast }$ is countably infinite. The length of the word u is denoted by $\left|u\right|$. By definition $\left|e\right|=0$. Obviously $\left|u\right|\in \mathbb{N}$ for each $u\ne e$, $\mathbb{N}=\{1,2,\dots \}$ stands for the set of natural numbers.

For $u\in X^{\ast }$ and $v\in Y^{\ast }$, if $\left|u\right|=\left|v\right|$ we write $\left(u\right|v)\in {\left(X\right|Y)}^{\ast }$ to distinguish it from the case $(u,v)\in X^{\ast }\times Y^{\ast }$. We denote by ${\left(X\right|Y)}^{\ast }$ the set of all input-output pairs of words of the same length:

$${\left(X\right|Y)}^{\ast }=\left\{\left(u\right|v)u\in X^{\ast },v\in Y^{\ast },\left|u\right|=\left|v\right|\right\}.$$

Definition 7.

Let $\mathcal{A}=(X,Q,Y,\mathcal{M})$ be FFMach.

For any $\left(u\right|v)\in {\left(X\right|Y)}^{\ast }$ the extended input-output behavior of $\mathcal{A}$ upon the law of composition * is determined by the square matrix $M\left(u\right|v)$ of order $\left|Q\right|$:

$$M\left(u\left|v\right.\right)=\left\{\begin{array}{c}M\left(x_1\left|y_1\right.\right)\ast \dots \ast M\left(x_k\left|y_k\right.\right),\hfill \\ if\left(u\right|v)=(x_1\dots x_k|y_1\dots y_k),k\ge 1\hfill \\ \\ E,if\left(u\right|v)=(e\left|e\right)\hfill \end{array}\right.,$$

(16)

where $M(x_1\left|y_1)\right.\ast \dots \ast M(x_k\left|y_k)\right.$, for $\ast =max-min$, $min-max$ or $max-$ products and $E=\left(e_{ij}\right)$ is the square matrix of order $\left|Q\right|$ with elements $e_{ij}$ determined as follows:

• if the composition is $max-min$ or $max-$ product,

  $$e_{ij}=\left\{\begin{array}{cc}1,\hfill & ifi=j\hfill \\ 0,\hfill & ifi\ne j\hfill \end{array}\right..$$
• if the composition is $min-max$,

  $$e_{ij}=\left\{\begin{array}{cc}0,\hfill & ifi=j\hfill \\ 1,\hfill & ifi\ne j\hfill \end{array}\right..$$

Each element $m_{q{q}^{\prime }}\left(u\right|v)$ in $M\left(u\right|v)$ is the degree of membership that $\mathcal{A}$ will enter state $q^{\prime }\in Q$ and produce output word $v\in Y^{\ast }$ under the input word $u\in X^{\ast }$ beginning at state $q\in Q$, after $\left|u\right|=\left|v\right|$ consecutive steps.

Definition 7 describes the extended input-output behavior of $max-min$, $min-max$ and $max-$product FFMachs.

4.2.2. Complete input-output behavior

When we consider FFMach as a ’black box’ we are not interested in the next state $q^{\prime }$. This is the essence of the input-output behavior of $\mathcal{A}$, determined by column-matrices $T{\left(u\right|v)}_{\left|Q\right|\times 1}$ as follows:

$$T{\left(u\right|v)}_{\left|Q\right|\times 1}=\left(t_q\left(u\right|v)\right)=\left\{\begin{array}{c}M\left(u\right|v)\ast E,if(u\left|v\right)\ne \left(e\right|e);\hfill \\ E,if\left(u\right|v)=(e\left|e\right),\hfill \end{array}\right.$$

(17)

where E is the $\left|Q\right|\times 1$ column-matrix with all elements equal to 1 if the composition is $max-min$ or $max-$product and with all elements equal to 0 if the composition is $min-max$.

Each element $t_q\left(u\right|v)$ of $T\left(u\right|v)$ in (17) determines the operation of $\mathcal{A}$ under the input word u beginning at state q and producing the output word v after $\left|u\right|=\left|v\right|$ consecutive steps.

For instance, if the FFMach is $max-min$ or $max-$product, the element

$$t_q\left(u\right|v)=\underset{q^{\prime }\in Q}{\vee }\left(m_{q{q}^{\prime }}\left(u\right|v)\right)$$

gives the way of achieving maximal degree of membership under the input word u, beginning at state q and producing the output word v.

We denote by $T_{\mathcal{A}}$ the complete behavior matrix of $\mathcal{A}$. It is semi-infinite matrix with $n=\left|Q\right|$ rows and with columns $T\left(u\right|v),$$\left(u\right|v)\in {\left(X\right|Y)}^{\ast }$, computed by (17) and ordered according to the lexicographical order in ${\left(X\right|Y)}^{\ast }$, see Table 2.

Let us mention that $t_q\left(e\right|e)=1$ for $i=1,...,n$ if FFMach is $max-min$ or $max-$product, $t_q\left(e\right|e)=0$ for $i=1,...,n$ if FFMach is $min-max.$

4.3. Behavior Matrix

For any FFMach $\mathcal{A}$ its complete input-output behavior matrix $T_{\mathcal{A}}$ is semi-infinite – it has finite number of rows (equal to the number of states in Q) and infinite number of columns. Since $T_{\mathcal{A}}$ is semi-infinite, one can not solve traditional problems – equivalence of states, reduction of states, minimization of states, because all of them explore $T_{\mathcal{A}}$.

The matrix $B_{\mathcal{A}}$ obtained by omitting all columns from $T_{\mathcal{A}}$ that are linear combination of the previous columns is called behavior matrix of $\mathcal{A}$.

We may compute the behaviour matrix for FFMach, but the time complexity of the algorithm is exponential. Algorithm 15 in subSection 5.8 extracts the finite behaviour matrix $B_{\mathcal{A}}$ from the complete behavior matrix $T_{\mathcal{A}}$, it implements linear dependency/independency from subSection 4.1.

4.4. Equivalence of States, Reduction, Minimization

$B_{\mathcal{A}}$ captures all the properties of $T_{\mathcal{A}}$ and provides solving equivalence, reduction and minimization problems for FFMachs.

For the FFMach $\mathcal{A}=(X,Q,Y,\mathcal{M})$ the states $q_i\in Q$ and $q_j\in Q$ are called equivalent if the input-output behavior of $\mathcal{A}$ when beginning with initial state $q_i$ is the same as its input-output behavior when beginning with initial state $q_j.$ It means that the $i-$th and $j-$th rows are identical in $T_{\mathcal{A}}$ and in $B_{\mathcal{A}}$ . Since $T_{\mathcal{A}}$ is semi-infinite, we can not derive equivalence of states from it. In order to solve this problem we first have to extract $B_{\mathcal{A}}$ from $T_{\mathcal{A}}$ and then to find identical rows in $B_{\mathcal{A}}$. Then we left only one of those identical rows.

FFMach $\mathcal{A}=(X,Q,Y,\mathcal{M})$ is in reduced form if there does not exist equivalent states in Q.

A FFMach $\mathcal{A}=(X,Q,Y,\mathcal{M})$ is not in minimal form if the input-output behavior of $\mathcal{A}$ when it begins with initial state $q_i$ is the same as its input-output behavior when it begins with initial distribution over Q, isolating the state $q_i.$ Formally it means that the $i-$th row of $T_{\mathcal{A}}$ ($B_{\mathcal{A}}$, respectively) is a linear combination of the other rows.

Software from [38] implements this idea: if the $i-$th row of $B_{\mathcal{A}}$ is a linear combination of the other rows of $B_{\mathcal{A}}$, then the FFMach $\mathcal{A}$ is not in minimal form.

Algorithms and software for obtaining the behaviour matrix, as well as for solving reduction and minimization problems are given in subSection 5.8 and Section 6.4 with working examples in Appendix A.3.

4.5. Linear Optimization

We present optimization problem when the linear objective function is

$$Z=\sum _{j=1}^n{c}_j{x}_j,c_j\in \mathbb{R},0\le x_j\le 1,1\le j\le n,$$

(18)

with traditional addition and multiplication; $c=(c_1,...,c_n)$ is the weight (cost) vector and the function (18) is subject to FLSEs or FLSIs as constraint

$$A{\ast }_{BL}X=B,orA{\ast }_{BL}X\le B,orA{\ast }_{BL}X\ge B.$$

(19)

The aim is to minimize or maximize Z subject to one of the constraints in (19).

The linear optimization problem was investigated first for $max-min$ composition by Fang and Li [49], Guu and Wu [50], Peeva [41]; for $max-$product composition by Loetamphong and Fang [51], Peeva and co-authors [52]; for $BL-$algebras by Peeva and Petrov [53]. If the objective function becomes $Z={max}_{j=1}^{n}min(c_j,x_j)$ with $c_j\in [0,1]$ the so called latticized linear programming problem was studied by Wang et al  [54].

Since the solution set of (19) is non-convex, traditional linear programming methods cannot be applied. In order to find the optimal solution for (18), the problem is converted into a $0-1$ integer programming problem. Fang and Li in  [49] solved this $0-1$ integer programming problem by branch and bound method with jump-tracking technique. Improvements of this method are proposed in [55] by providing fixed initial upper bound for the branch and bound part and with update of this initial upper bound. Conventional investigations  [49],  [50],  [51],  [55] begin with various attempts to find some feasible solutions, then to find an arbitrary solution of the optimization problem subject to constraint from (19) and proceed with it in the solution set for search of minimal solutions. The procedures are clumsy, the main obstacle is how to solve (19).

The algorithms and software from [38] proceeds for solving optimization problems, among them are:

i)

Minimize the linear objective function (18), subject to FLSEs constraints, when the composition is $max-min$, $max-$product or $max-$ukasiewicz;

ii)

Maximize the linear objective function (18), subject to FLSEs, when the composition is $max-min$, $max-$product or $max-$ukasiewicz;

iii)

Minimize or maximize the linear objective function (18), subject to FLSEs constraints, when the composition is $min-max.$

For clarity of exposition we present the ideas for $max-min$ and $max-$product case.

4.6. Minimize the Linear Objective Function, Subject to FLSEs Constraints, When the Composition is $max-min$ or $max-$Product Case.

We present how to minimize the linear objective function (18), subject to constraints: a fuzzy linear system of equations with $max-min$ or $max-$ product composition

$$A\circ X=B,$$

(20)

where $A={\left(a_{ij}\right)}_{m\times n}$ with $a_{ij}\in [0,1]$ stands for the matrix of coefficients, $X={\left(x_i\right)}_{n\times 1}$, $x_i\in [0,1]$ stands for the matrix of unknowns, $B={\left(b_i\right)}_{m\times 1}$, $b_i\in [0,1]$ is the right-hand side of the system, the $max-min$ or $max-$ product composition is written as ∘ and $c=(c_1,...,c_m)$ is the weight (cost) vector.

We apply the software from  [38] for computing the greatest solution and all minimal solutions of the consistent system $A\circ X=B$. Then we solve the linear optimization problem, first decomposing Z into two vectors $Z^{\prime }=(c_1^{\prime },c_2^{\prime },...,c_n^{\prime })$ and $Z^{\prime \prime }=(c_1^{\prime \prime },c_2^{\prime \prime },...,c_n^{\prime \prime })$, such that $Z=Z^{\prime }+Z^{\prime \prime }$ and:

$$c_i=c_i^{\prime }+c_i^{\prime \prime },foreachi=1,...,n,$$

$$c_i^{\prime }=\left\{\begin{array}{c}{c}_i,if{c}_i\ge 0,\\ 0,if{c}_i<0\end{array}\right.,$$

(21)

$$c_i^{\prime \prime }=\left\{\begin{array}{c}0,if{c}_i\ge 0,\\ c_i,if{c}_i<0\end{array}\right..$$

(22)

Now the original problem: to minimize Z subject of constraints (20), is split into two problems, namely to minimize both

$$Z^{\prime }=\sum _{i=1}^m{c}_i^{\prime }{x}_i$$

(23)

and

$$Z^{\prime \prime }=\sum _{i=1}^m{c}_i^{\prime \prime }{x}_i$$

(24)

with constraints (20).

Z takes its minimum when $Z^{\prime }$ takes its minimum and $Z^{\prime \prime }$ takes its maximum. Hence for the problem (24) the optimal solution is $\widehat{X}=(\widehat{x_1},...,\widehat{x_n})=X_{gr}$, and for the problem (23) the optimal solution $\breve{X}=(\breve{x_1},...,\breve{x_n})$ is among the minimal solutions of the system (20). In this case the optimal solution of the problem is $X^{\ast }=(x_1^{\ast },...,x_n^{\ast })$, where

$$x_i^{\ast }=\left\{\begin{array}{c}\widehat{x_i},if{c}_i<0\\ \breve{x_i},if{c}_i\ge 0\end{array}\right..$$

(25)

4.7. Maximize the Linear Objective Function, Subject to Constraints $A\circ X=B$, When the Composition is $max-min$ or $max-$Product

In order to maximize the linear objective function Z, we again split it, but now for (24) the optimal solution is among the minimal solutions of the system (20), for (23) the optimal solution is $X_{gr}$. In this case the optimal solution is $X^{\ast }=(x_1^{\ast },...,x_n^{\ast })$, where

$$x_i^{\ast }=\left\{\begin{array}{c}\breve{x_i},if{c}_i<0\\ \widehat{x_i},if{c}_i\ge 0\end{array}\right..$$

(26)

4.8. Algorithm for Finding Optimal Solutions

We propose the algorithm for finding optimal solutions if the linear objective function is Z and the composition is $max-min$ or $max-$product for the constraints (20).

In any of these cases the optimal value is

$$Z^{\ast }=\sum _{i=1}^n{c}_i{x}_i^{\ast }$$

(27)

The algorithm for finding optimal solutions is based on the above results.

Algorithm 4 Algorithm for finding optimal solutions

Step 1: Enter the matrices $A_{m\times n}$, $B_{m\times 1}$ and the cost vector $C_{1\times n}$. Step 2: Implement the software from [38] to establish consistency of the system (20). If the system is inconsistent - Exit. Step 3: Implement the software from [38] to compute $X_{\mathrm{gr}}$ and all minimal solutions of (20). Step 4: If finding $Z_{max}$ go to Step 6. Step 5: For finding $Z_{min}$ compute $x_i^{\ast }$, $i=1,...,n$ according to (25). Go to Step 7. Step 6: For finding $Z_{max}$ compute $x_i^{\ast }$, $i=1,...,n$ according to (26). Step 7: Compute the optimal value according to (27).

Algorithms and software for the above described linear optimization problems are given in subSection 5.9 and Section 6.5 with working examples in Appendix A.4.

5. Algorithms

In this section, we present algorithms that have been developed over the years and are used in the software [38]. These algorithms cover fuzzy matrix compositions, inverse problem resolution, linear dependence and independence, fuzzy optimization, and finite fuzzy machines. They are designed to support a variety of fuzzy algebras, including max-min, min-max, max-product, and all the others defined in Definitions 1 - 3.

However, in order to keep the exposition focused and accessible, we restrict attention in this chapter to the max-min case. This is justified not only by its historical role in the development of fuzzy relational methods, but also by the fact that it forms the basis of many practical applications. Other compositions introduced in Definitions 1 - 3 will be demonstrated in the following Chapter Section 6.

The inverse problem for fuzzy linear systems is particularly challenging when computing the complete solution set. The algorithms presented here incorporate structural optimizations to address this challenge while preserving the mathematical properties of the underlying systems.

5.1. Fuzzy Matrix Compositions

Direct problem in FRC means fuzzy matrix composition. According to Definitions 1 - 3, fuzzy matrix composition is defined by applying two fuzzy operations. Similarly to conventional matrix multiplication, the general procedure computes each entry $c_{ij}$ of the resulting matrix C from the row i of matrix A and column j of matrix B.

Algorithm 5 Generic Fuzzy Matrix Composition

Step 1: Input fuzzy matrices A, B, and the pair of fuzzy operations Step 2: For each entry $c_{ij}$ in the resulting matrix C: Step 3: For each index k: Step 4: Apply the first operation to $a_{ik}$, $b_{kj}$ Step 5: Aggregate all intermediate results using the second operation Step 6: Output matrix C as the composition between A and B

Complexity: Algorithm 5 computes each entry $c_{ij}$ by evaluating k pairs and aggregating the result, leading to a total of $mnk$ basic operations. This gives a complexity of $O\left(mnk\right)$.

This algorithm serves as a general template for all supported compositions. Each specific case, such as max-min, min-max, max-product, or compositions based on fuzzy implications is implemented as an alias function that calls this general routine with the corresponding pair of fuzzy operations.

5.2. Inverse Problem Resolution for Fuzzy Linear Systems of Equations

The inverse problem in FRC consists of solving an FLSEs of the form $AX=B$ or $A{\ast }_{BL}X=B$ or $A{\to }_{BL}X=B$ where A and B are known fuzzy matrices, X is the unknown, and the composition is defined as in Definitions 1 - 3. Solving this problem means finding the complete solution set $\mathbb{X}$ under the selected fuzzy composition.

While direct composition is computationally straightforward, solving the inverse problem, especially for obtaining the complete solution set, is significantly more complex. In what follows, we focus on the classical case of max-min composition, where most theoretical and algorithmic developments are concentrated.

Over the years, several different algorithms have been developed and implemented for solving the inverse problem in Fuzzy Calculus Core. These earlier methods, based on the algebraic-logical framework and analytical constructions, are described in detail in Section 3. The current version of the software incorporates an optimized algorithm that builds on these foundational techniques.

The method described below represents the latest and most efficient version currently implemented in the software. It is designed to compute the greatest solution and recursively generate all minimal solutions through domination and recursive search.

Algorithm Overview

The algorithm for solving FLSEs $A•X=B$ under max-min composition is divided into the following stages:

• Stage 1: Initializations
• Stage 2: Find the greatest solution
• Stage 3: Domination
• Stage 4: Recursive extraction of all lower solutions

Stage 1: Initializations

Algorithm 6 Stage 1: Initializations

Step 1: Enter fuzzy matrix A and fuzzy vector B Step 2: Initialize help matrix H with the same size as A, filled with 0 Step 3: for each $a_{ij}$ in A do Step 4:     if $a_{ij}\ge b_i$ then Step 5:         Set $h_{ij}=b_i$ Step 6:     end if Step 7: end for

Complexity: Algorithm 6 has complexity $O\left(mn\right)$, where m and n are the number of rows and columns, respectively, of the input matrix A.

Stage 2: Find the greatest solution

Earlier algorithms for computing the greatest solution relied on slow residuum-based expressions or complicated algorithms based on S, E, G, and H elements, presented in Section 3.4.

The current implementation avoids implication operators entirely. Instead, it computes the greatest solution directly from the help matrix using column operations. The consistency of the system is checked during the same pass, marking all covered rows without additional matrix traversal. This leads to a significant improvement in performance and simplicity.

Algorithm 7 Stage 2: Find the greatest solution

Step 1: Initialize vector X with ones representing the greatest solution Step 2: Initialize vector $IND$ with zeros to mark system consistency Step 3: for each j in X do Step 4:     Set $x_j$ to the minimum of all $h_{ij}$ such that $h_{ij}=b_i$ and $a_{ij}\ne b_i$ Step 5:     Mark all $IN{D}_i$ with 1 where $h_{ij}=b_i$ Step 6:     Set $h_{ij}=0$ for all i where $h_{ij}>x_j$ Step 7: end for Step 8: if there exists any 0 in $IND$ then Step 9:     System is inconsistent. Step 10: else Step 11:     System is consistent, its greatest solution is X Step 12: end if

Complexity: Algorithm 7 requires at most $mn$ comparisons in its worst-case scenario, giving $O\left(mn\right)$ complexity. This algorithm has better scenarios, where the solution can be found early for each variable, resulting in the best-case performance of $O\left(n\right)$.

Stage 3: Domination

Definition 8.

[15] Let $h_l=\left(h_{lj}\right)$ and $h_k=\left(h_{kj}\right)$ be the $l-$th and the $k-$th rows, respectively, in the help matrix H. If for each j, $1\le j\le n$, $h_{lj}\le h_{kj}$, then:

• $h_l$ is said to be a dominant row to $h_k$ in H
• $h_k$ is redundant row and can be removed from H

In earlier implementations, dominance was analyzed using a separate dominance matrix. In the current version, this step is evaluated directly over the help matrix.

Algorithm 8 Stage 3: Domination

Step 1: for each row i in H do Step 2:     for each row $j\ne i$ in H do Step 3:         if row j dominates row i then Step 4:            Mark row i for deletion Step 5:            Break Step 6:         else if row i dominates row j then Step 7:            Mark row j for deletion Step 8:         end if Step 9:     end for Step 10: end for Step 11: Remove all marked rows from H

Complexity: Algorithm 8 is giving a worst-case complexity of $O\left(m^{2}n\right)$. In favorable cases where domination is detected early, many comparisons are skipped via early termination, resulting in best-case performance of $O\left(mn\right)$.

Stage 4: Recursive extraction of all lower solutions

The current implementation uses a recursive procedure that incrementally builds each solution from valid components in the help matrix. At every step, the algorithm maintains a set of uncovered rows and explores only those elements that contribute to the current potential solution. Branches that cannot lead to a complete solution are skipped early. This method avoids full enumeration of all combinations, and provides a conceptually simpler and more memory-efficient procedure.

Algorithm 9 Stage 4: Recursive extraction of all lower solutions

Step 1: Initialize empty set of solutions $\mathbb{X}$ Step 2: Initialize candidate solution X with zeros Step 3: Initialize marking vector M with zeros, meaning no rows are initially marked Step 4: Start recursion with the first unmarked row $i=1$ Step 5: for each column j such that $h_{ij}>0$ do Step 6:     Set $x_j=h_{ij}$ Step 7:     Mark all rows k such that $h_{kj}>0$ by setting $m_k=1$ Step 8:     if all rows are marked then Step 9:         if X is not dominated by any other $Y\in \mathbb{X}$ then Step 10:            Remove all $Y\in \mathbb{X}$ such that $X\le Y$ Step 11:            Add X to $\mathbb{X}$ Step 12:         end if Step 13:     else Step 14:         Let $i^{\prime }$ be the next unmarked row in M Step 15:         Repeat this loop recursively with updated X, M, and $i^{\prime }$ Step 16:     end if Step 17: end for

Figure 1. Visual representation of Algorithm 9.

Figure 1. Visual representation of Algorithm 9.

Complexity: The time complexity of the recursive extraction of all lower solutions from H is exponential, but in this case, it depends on the number of non-zero elements in the help matrix, which itself reflects the number of possible lower solutions. More precisely, in the worst case, the number of recursive paths can grow to $O\left(d^m\right)$, where d is the average number of non-zero entries per row, and m is the number of rows in H which may be smaller than in A, due to the domination step. In practice, however, the algorithm is highly efficient due to early pruning of incomplete branches and the reduced size of H after domination.

5.3. Complete Algorithm for Solving $A•X=B$

The previous procedures are combined into a complete algorithm for solving FLSEs $A•X=B$.

Algorithm 10 Solving the inverse problem for $A•X=B$

Step 1: Initialize the help matrix H using Algorithm 6 Step 2: Compute the greatest solution and check consistency using Algorithm 7 Step 3: if the system is consistent then Step 4:     Perform row domination on H using Algorithm 8 Step 5:     Extract all lower solutions using Algorithm 9 Step 6: end if

Complexity: The total complexity is dominated by the recursive enumeration of all lower solutions (Algorithm 9), which is exponential in the number of rows in the help matrix. All preceding steps—help matrix construction, consistency check, and row domination—run in $O\left(mn\right)$ or $O\left(m^{2}n\right)$ time.

The presented algorithm can be extended to other types of systems and compositions. Similar procedures are implemented in [38] for all composition operators defined in Definition 1, illustrating the potential for both algorithmic and theoretical generalizations of this approach.

5.4. Inverse Problem Resolution for Fuzzy Linear Systems of ≥ Inequalities

While FLSEs of the form $A•X=B$ have been widely studied, FLSIs such as $A•X\ge B$ have received comparatively little attention in the literature. To the best of our knowledge, no other general-purpose software tools for solving such systems exist.

The algorithms presented in the previous section for solving FLSEs can be extended to handle systems of inequalities with only minor modifications. In this section, we focus on the case $A•X\ge B$, where the greatest solution is trivial if it exists, but lower solutions must still be extracted recursively. We cannot assume that the corresponding system of equations has a solution, there are cases where the inequality system is consistent, even though the equation system is not. We highlight the key differences in structure and interpretation compared to the case of equations.

Algorithm 11 Solving the inverse problem for $A•X\ge B$

Step 1: Use Algorithm 6 to find the help matrix H Step 2: Initialize vector $IND$ with zeros Step 3: Mark all $IN{D}_i$ with 1 where $h_{ij}=b_i$ Step 4: if any element of $IND$ is still zero then Step 5:     The system is inconsistent Step 6: else Step 7:     Set X with ones for the greatest solution Step 8:     Apply Algorithm 8 to eliminate dominated rows from H Step 9:     Use Algorithm 9 to compute all lower solutions Step 10: end if

Complexity: The computational complexity is identical to that of the algorithms described in Section 5.2.

5.5. Inverse Problem Resolution for Fuzzy Linear Systems of ≤ Inequalities

FLSIs of the form $A•X\le B$ are structurally simpler than equations or ≥ inequalities. A trivial lower solution X of all zeros always exists, and no recursive extraction is required. The problem reduces to finding the greatest solution that satisfies all inequality constraints. This is performed using the same method as in the equality case, with minor simplifications.

Algorithm 12 Solving the inverse problem for $A•X\le B$

Step 1: Use Algorithm 6 to find the help matrix H Step 2: Initialize vector $IND$ with zeros Step 3: Mark all $IN{D}_i$ with 1 where $h_{ij}=b_i$ Step 4: if any element of $IND$ is still zero then Step 5:     The system is inconsistent Step 6: else Step 7:     Use Algorithm 7 to find the greatest solution X Step 8:     Set X zeros as the trivial lower solution Step 9: end if

Complexity: The complexity is significantly lower than in Section 6.2. No domination or recursion is required. We are left with just the greatest solution, which is computed in $O\left(mn\right)$ time.

5.6. Fuzzy Linear Combination of Vectors

A fuzzy vector B is said to be a linear combination of the columns of a matrix A under a fixed composition (e.g. max-min) if there exists a fuzzy vector X such that $A•X=B$. This is equivalent to solving the inverse problem.

To determine whether such a combination exists, we reuse the algorithms for computing the greatest solution described in Section 5.2. Since we are only interested in the existence of a solution, it is sufficient to validate that the system is consistent i.e. there is no need to extract all lower solutions.

Algorithm 13 Check if B is a fuzzy linear combination of the columns of A

Step 1: Enter fuzzy matrix A and vector B Step 2: Check if the system $A•X=B$ is consistent according to Algorithm 7 Step 3: if the system is consistent then Step 4:     B is a fuzzy linear combination of the columns in A Step 5: else Step 6:     B is not a fuzzy linear combination of the columns in A Step 7: end if

Complexity: The complexity is the same as the consistency check in Algorithm 7.

5.7. Fuzzy Linear Independence

A set of fuzzy vectors $\{A_1,A_2,\dots ,A_n\}$ is said to be linearly independent if no vector in the set can be expressed as a fuzzy linear combination of the others vectors under a fixed composition (e.g. max-min).

Our implementation tests linear independence by iteratively checking, each column $A_i$ in a fuzzy matrix A, for a linear combination of the remaining columns. This is done using the procedure from Algorithm 13.

Algorithm 14 Check if the columns of A are fuzzy linearly independent

Step 1: Enter fuzzy matrix A Step 2: Initialize empty list D to hold dependent columns Step 3: for each column $A_i$ in A do Step 4:     Let $A_i^{\ast }$ be the matrix A without column $A_i$ Step 5:     Check if $A_i$ is a fuzzy linear combination of $A_i^{\ast }$ using Algorithm 13 Step 6:     if $A_i$ is a linear combination of the other columns then Step 7:         Add index i to list D Step 8:     end if Step 9: end for Step 10: ifD is empty then Step 11:     The columns of A are fuzzy linearly independent Step 12: else Step 13:     D holds the list of all dependent columns Step 14: end if

Complexity: The algorithm performs n consistency checks, each with complexity $O\left(mn\right)$. Therefore, the total complexity is $O\left(m{n}^2\right)$.

5.8. Finite Fuzzy Machines

This subsection focuses on computing the behavior of max-min FFMach, see subSection 4.2. The behavior matrix is obtained by composing all reachable transitions iteratively until stabilization. The process terminates automatically when further compositions no longer change the result.

Here we present the algorithms used in [38] for computing the behavior matrix, and for two types of postprocessing: reduction and minimization presented in Section 4.2.

Algorithm 15 Compute behavior matrix B of a FFMach

Step 1: Enter the initial set of fuzzy transition matrices $\mathcal{M}=\{M_1,\dots ,M_k\}$ Step 2: Initialize the behavior matrix B with a single column of ones Step 3: repeat Step 4:     Let n be the current number of columns in B Step 5:     for $j=1$ to n do Step 6:         for each $M_i\in \mathcal{M}$ do Step 7:            Compute $b=M_i•B(:,j)$ Step 8:            if b is not a linear combination of the columns in B then Step 9:                Append b as a new column to B Step 10:            end if Step 11:         end for Step 12:     end for Step 13: until no new columns are added to B

Complexity: Let N denote the final number of behavioral vectors in the matrix B. At each iteration, up to $kn$ new candidates are generated, where k is the number of fuzzy transition matrices and n is the current number of columns in B. Each new vector is tested for fuzzy linear combination against the existing ones, using Algorithm 13, which has complexity $O\left(mn\right)$ for vectors of length m. Assuming no vector is rejected as dependent, the total number of candidates grows quadratically and each test becomes increasingly more expensive. The resulting worst-case complexity is approximately $O\left(km{N}^3\right)$.

Algorithm 16 Reduce the FFMach

Step 1: Enter the behavior matrix B Step 2: Identify all duplicate columns in the transposed matrix $B^T$ Step 3: Remove the corresponding rows from the original matrix B

Complexity: The software [38] uses MATLAB’s built-in unique function, which is based on lexicographical sorting. The overall complexity is approximately $O(mlogm)$, where m is the number of rows.

Algorithm 17 Minimize states of the FFMach

Step 1: Enter the behavior matrix B Step 2: Use Algorithm 14 to identify all dependent columns in the transposed matrix $B^t$ Step 3: Remove the corresponding rows from the original matrix B

Complexity: The algorithm internally applies Algorithm 14, and therefore has the same complexity.

5.9. Fuzzy Optimization

Let us remind (see Section 4.5) that the linear objective function

$$Z=\sum _{j=1}^n{c}_j{x}_j$$

is subject to fuzzy constraints of the form

$$A•X=B,A•X\le B,\mathrm{or}A•X\ge B.$$

The algorithm splits the objective function into two parts: $Z^{\prime }$, corresponding to components with $c_i\ge 0$, and $Z^{\prime \prime }$, corresponding to components with $c_i<0$. When minimizing Z, the variables in $Z^{\prime }$ are selected among all lower solutions, while the variables in $Z^{\prime \prime }$ are taken from the greatest solution of the constraint system. When maximizing - variables in $Z^{\prime }$ are taken from the greatest solution, and those in $Z^{\prime \prime }$ are selected among the lower solutions. The result is a combined vector X formed by selecting components from the appropriate solutions.

Since any value between a lower and an upper solution is itself a solution of the fuzzy system (see Section 3.2), the obtained X is also a valid solution of the system.

Algorithm 18 Fuzzy optimization under max–min FLSE constraints

Step 1: Enter fuzzy system $A•X★B$, where $★\in \{=,\le ,\ge \}$ Step 2: Enter objective weights C and optimization goal: minimize or maximize Step 3: Solve the constraint system using Algorithm 10, 12, or 11, depending on ★ Step 4: Split the objective function into $Z^{\prime }$ (for $c_j\ge 0$) and $Z^{\prime \prime }$ (for $c_j<0$) Step 5: if goal is minimize then Step 6:     Take the corresponding components for $Z^{\prime \prime }$ from the greatest solution Step 7:     Evaluate $Z=Z^{\prime }+Z^{\prime \prime }$ for all lower solutions Step 8:     For $Z^{\prime }$, select the solution that minimizes Z Step 9: else  if goal is maximize Step 10:     Take the corresponding components for $Z^{\prime }$ from the greatest solution Step 11:     Evaluate $Z=Z^{\prime }+Z^{\prime \prime }$ for all lower solutions Step 12:     For $Z^{\prime \prime }$, select the solution that maximizes Z Step 13: end if

Complexity: The algorithm internally applies Algorithm 10, 12, or 11, depending on the constraint type, and inherits their respective complexity.

6. Software

This chapter presents the software implementation of the Fuzzy Calculus Core package [38]. The software provides tools for working with fuzzy matrices, solving fuzzy relational equations and inequalities, performing fuzzy optimization, and reduction and minimization of finite fuzzy machines.

The architecture of the system reflects the algorithmic structure introduced in Chapter Section 5. Each algorithm described there corresponds to a specific implementation in the respective module. Examples illustrating the use of each functionality are provided in Appendix A.

6.1. Architecture

The software is structured as a modular MATLAB package with a clear separation between low-level data representation and high-level problem-solving routines. Its architecture reflects the logical layering of the algorithms described in Chapter Section 5, and supports reuse, extension, and composition.

At the core of the system lies the fuzzyMatrix module, which implements data structures and operations for fuzzy matrices and vectors. On top of this base layer are modules for solving fuzzy linear systems of equations and inequalities (fuzzySystem), for solving optimization problems (fuzzyOptimizationProblem), and for finding behavior, minimizing and reducing fuzzy finite machines (fuzzyMachine).

An overview of the module dependencies is presented in Figure 2.

The system supports a variety of fuzzy matrix compositions, as defined in Definitions 1–3. These are:

• Max-min
• Min-max
• Max-product
• Max-$ϵ$
• Łukasiewicz composition
• Min-${\to }_G$ (Gödel implication)
• Min-${\to }_P$ (Product implication)
• Min-${\to }_L$ (Łukasiewicz implication)

All compositions are supported uniformly across the system: they can be used not only for direct matrix operations, but also in solving both the direct and the inverse problems, as well as in fuzzy optimization and fuzzy finite machines.

6.2. The fuzzyMatrix Module

The fuzzyMatrix module provides the core data structure for representing fuzzy matrices. It is implemented as a MATLAB class that extends the built-in double type, enabling native support for arithmetic operations while enforcing fuzzy semantics.

A fuzzyMatrix object can be initialized in one of the following ways:

• Without arguments, creating an empty fuzzy matrix.
• With a single argument, which must be a real matrix with values in the $[0,1]$ interval.
• With two scalar arguments, interpreted as dimensions, which creates a zero-filled matrix of the corresponding size.

Internally, all matrix elements are stored as standard numerical arrays. The class overrides basic operations such as addition and subtraction to ensure that results remain within the fuzzy domain.

In addition to standard elementwise operations and transposition, the module implements a wide range of fuzzy matrix compositions:

• maxmin, minmax, maxprod
• minalpha, maxepsilon, mindiamond
• godel, goguen, lukasiewicz, maxlukasiewicz
• minprobabilistic, minbounded, maxdelta, maxgama

Each of these is implemented through the general-purpose method fcompose, which composes two fuzzy matrices using a pair of operations.

The module also provides methods for fuzzy linear combination (is_lincomb) and fuzzy linear independence (is_linindep), allowing users to test whether a given vector can be expressed as a fuzzy combination of others under a specified composition.

Examples illustrating the use of this module can be found in Appendix A.1.

6.3. The fuzzySystem Module

The fuzzySystem module provides a high-level abstraction for solving fuzzy relational problems of the form =, ≤, and ≥. Both the direct and the inverse problems are supported, under all compositions listed in Section 6.1.

The module acts as a front-end interface to the lower-level algorithms described in Chapter Section 5. A fuzzySystem object is constructed by specifying the composition type and the equality/inequality type, together with two fuzzyMatrix instances. It delegates the computation to the appropriate method depending on the composition.

For direct problems (computing B), the module internally calls the fcompose method from the fuzzyMatrix class. For inverse problems the underlying logic follows Algorithms  10,  11 and  12, depending on the selected system type.

The main interface function solve_inverse returns a structure with standardized fields, which include:

• gr – a set of upper solutions, which contains either the greatest solution or all upper solutions, depending on the composition
• low – a set of lower solutions, which contains either the lowest solution or all lower solutions;
• exist – a boolean flag indicating whether the system is compatible or incompatible.

In certain situations, it is sufficient to check whether a system is consistent or to obtain a single representative solution. In such cases, the algorithm computes only the greatest (resp. lowest) solution, avoiding the extraction of the full solution set. This mode is significantly more efficient and is based on Algorithm 7, which also serves as a consistency check. When full enumeration is needed, setting full=true triggers the recursive Algorithm 10 form Section 5.

All input matrices are instances of fuzzyMatrix, and all computed results are returned in the same class.

Examples of usage are provided in Appendix A.2.

6.4. The fuzzyMachine Module

The fuzzyMachine module implements operations for finding behavior, reduction and minimization of FFMach. A FFMach is defined by a set of fuzzy matrices, one for each input/output pair of words. These matrices are stored in a cell array and passed to the constructor. The user specifies the composition type, the maximum word length to consider, and whether to return the full behavior matrix, a reduced or a minimized one.

The main method find_behavior() applies the matrices to compute the machine’s behavior. If postprocessing is enabled, the resulting behavior matrix is reduced or minimized by eliminating redundant columns or states. The minimization procedure removes states whose behavior vectors can be expressed as fuzzy linear combinations of others.

All compositions listed in Section 6.1 are supported. The final behavior matrix is stored internally and can be accessed or exported for further analysis.

Examples are provided in Appendix A.3.

6.5. The fuzzyOptimizationProblem Module

The fuzzyOptimizationProblem class provides a unified interface for solving fuzzy linear optimization problems of the form described in Section 4.5. A problem instance is created by specifying:

• an objective weight (cost) vector C;
• a system of constraints defined by a fuzzySystem object and supporting all of the fuzzySystem options such as compositions, and inequalities

Internally, the functions first solve the inverse problem for the fuzzySystem object to obtain the complete solutions set for the system of constraints. Then, it either evaluates the objective function on all solutions, according on the optimization goal.

The methods minimize() and maximize() fill the object_solution and the object_value properties of the fuzzyOptimizationProblem object containing respectively, the optimal fuzzy vector x and the objective value z for the optimization problem. All properties of the system of constrains are held in the fuzzySystem object (e.g., consistency status or complete solutions set).

Examples of usage are provided in Appendix A.4.

7. Conclusions

In 1977 Peeva finished her PhD Thesis on categories of stochastic machines. The main attention was paid on computing behavior, establishing equivalence of states and equivalence of stochastic machines, as well as on reduction and minimization. All these problems were solved for stochastic machines, using linear algebra: they require solving linear system of equations with traditional algebraic operations, establishing linear dependence or linear independence of vectors, Noetherian property. After 1980 she was interested in similar problems, but for finite fuzzy machines. Nevertheless that stochastic machine seems to be similar to fuzzy machine, the main obstacle was that linear algebra and fuzzy algebras are completely unrelated. In order to investigate fuzzy machines, two types of algebras, called $max-min$ algebra and $max-$product algebra have been developed. The role played by these algebras in the theory of $max-min$ and $max-$product FFMachs is supposed to be the same as that played by linear algebra in the theory of stochastic machines. But these algebras propose tremendous manipulations for solving problems for FFMachs. Supplementary, there are several issues in [56], that are either incomplete or confused. This motivated the author to develop theory, algorithms and software for FFMachs. Obviously developing theory and software for solving problems for FFMachs requires first to develop theory and software for FRC and then to implement it for FFMachs. At this time these problems were open and she intended first to develop method, algorithm and software for solving FLSEs when the composition is $max-min$, then to implement it to FFMachs. In 1980 she had no idea whether the solution exists and how long it will take her to solve the problems. Just now, after 40 years, we have positive answers to all these questions.

She worked on this subject with great pleasure during the years with her bachelor, master and PhD students, now Prof. Dr. Yordan Kyosev, Assoc. Prof. Dr. Zlatko Zahariev, Dobromir Petrov and Ivailo Atanassov to whom she would like to express deep and heartfelt gratitude.

Zlatko Zahariev joined this work in 2005 as a PhD student under the supervision of Prof. DSc. Ketty Peeva. His main focus was on developing algorithmic versions the methods, creating their software implementation, and contributing to the ongoing effort toward designing unified algorithms for solving the problems described in this paper. He considers himself fortunate to have the opportunity to work under Prof. DSc. Peeva’s guidance, and is deeply grateful to her for the years of meaningful and fulfilling scientific collaboration.

The authors thank the journal team for the special offer to write this article.

Appendix A.1. Examples for the fuzzyMatrix module

This section provides usage examples for the fuzzyMatrix module, including initialization, basic operations, and several types of fuzzy matrix composition.

Appendix A.2. Examples for the fuzzySystem module

This section provides example uses of the fuzzySystem class, illustrating how to construct and solve fuzzy relational systems under various compositions and inequality types.

Appendix A.3. Examples for the fuzzyMachine module

Appendix A.4. Examples for the fuzzyOptimizationProblem module

This section provides examples demonstrating the use of the fuzzyOptimizationProblem class to define, minimize and maximize fuzzy optimization problems.