Likelihood - estimation on parameters I make confusions in the using of (1) Maximum Likelihood to find (approximately) the original signal knowing the data observed and (2)the using of Maximum Likelihood to find estimations of parameters of the PSF.
Goal : Find (up to some point) the original signal :
I start from this general definition (in discretized) : $$y=H\,x+w\quad(1)$$ (with $w$ a white noise)
Question: How can I demonstrate this relation ? (it seems that we should start from a discrete convolution product, doesn't it ? Then, the correct expression would be rather : $y=H*x+w$ with $*$ the product covolution
Thanks
 A: When the image is a point source.
The convolution of the PSF with some original image becomes a simple function when the original image is a point source (ie it can be modelled with a delta function).
$$\begin{align}
d_{\text{after}}(r,c) & =  (d_{\text{before}}* PSF)(r,c) \tag{1a}\label{eq1a} \\
&= \int \int d_{\text{before}}(x,y) PSF(r-x,c-y) dx dy \tag{1b}\label{eq1} \\
& =  \int \int \delta(r_0-x,c_0-y) PSF(r-x,c-y) dx dy \tag{1c}\label{eq1c} \\
& =  PSF(r-r_0,c-c_0) \tag{1d}\label{eq1d}
\end{align} $$
see for more information here: https://en.wikipedia.org/wiki/Dirac_delta_function#Translation and https://en.wikipedia.org/wiki/Convolution#Translational_equivariance
So, in the case of a star/point-source, then your image is just the (single) PSF translated to the position of the coordinates of the star. You can forget about the convolution (in the sense of intergrals) from this point and just work with the equation:
$$d(r,c) =  a \cdot PSF_{\alpha,\beta}(r-r_0,c-c_0) + b + \epsilon(r,c) \tag{2}\label{2}$$
Your deconvolution problem is like you have some observation $d(r,c)$ that is dependent on a star positioned at $r_0$, $c_0$, that gives a signal with amplitude $a$ and is spread out by some PSF that depends on $\alpha$, $\beta$, and where you have an additional (uniform?) background $b$ and white Gaussian noise $\epsilon(r,c)$ with variance $\sigma_\epsilon^2$.
When the other variables $\mathbf{r_0,c_0,\alpha,\beta}$ are known
When the other variables $\mathbf{r_0,c_0,\alpha,\beta}$ are known then the function $PSF_{\alpha,\beta}(r-r_0,c-c_0)$ is fixed for every pair $r,c$.
Then this becomes like a ordinary least squares model (OLS). $$ y_{r,c} = b+ a 
   x_{r,c}  + \epsilon_{r,c} \tag{3}\label{eq3}$$ with $$\begin{array}{rcl} y_{r,c} &=& d(r,c) \\ x_{r,c} &=& PSF_{\alpha,\beta}(r-r_0,c-c_0) \\ \epsilon_{r,c} &=& \epsilon(r,c) 
    \end{array}$$
When the other variables are not known
When the other variables are not known then it becomes like a non-linear regression problem, and you could solve it by a grid search or by using some algorithm for solving non-linear least squares prolems (of which there are many).

What is the deal with the matrix equation?
You seem to be looking for some equivalence between convolution $f*h$ and some matrix equation $f.h$. But you should forget about this convolution since we are dealing with the simpler equation (2) instead of (1a).
One thing that might resemble an equivalence might be the fact that a convolution is equivalent to a product in frequency space, and it is this fact that you might use when you de-convolving when you are not dealing with point sources (see for more information https://en.wikipedia.org/wiki/Deconvolution but note that this is not what you are doing in your exercise, which is a simpler case).
You can see the expression that you are working with (which would become a sum of PSFs for multiple point sources) as a sum of PSF with a fixed number of parameters. So you can put it in a form where you only have a limited number of parameters to estimate.
On the other hand, in the case of the convolution where the source image is not a point source (not limited numbers) the source image will be defined by an infinite sum of PSF functions and the system is underdetermined. You can not solve it as a simple regression problem where you minimize the sum of squares of the error.

Answering/helping some of your many questions
1

(it seems that we should start from a discrete convolution product,
doesn't it ? Then, the correct expression would be rather : $y=H∗x+w$
with $∗$ the product convolution )

In your particular question the star is considered a point source. In this case you do not need to use a convolution (or you do a convolution, but it is with a delta function see sifting or sampling property)
2

What is the link between H and the PSF used above (Moffat PSF) ?

In matrix form you have the function $d(r,c) =  a \cdot PSF_{\alpha,\beta}(r-r_0,c-c_0) + b + \epsilon(r,c)$ as:
$$\begin{bmatrix}d(1,1) \\ d(1,2) \\ d(1,3) \\ \vdots \\ d(20,20)  \end{bmatrix} = \begin{bmatrix} 1 & PSF_{\alpha,\beta}(1-r_0,1-c_0) \\ 1 & PSF_{\alpha,\beta}(1-r_0,2-c_0)  \\ 1 & PSF_{\alpha,\beta}(1-r_0,3-c_0)  \\ \vdots & \vdots \\ 1 & PSF_{\alpha,\beta}(20-r_0,20-c_0)  \end{bmatrix} \times  \begin{bmatrix}b \\ a \end{bmatrix} + \begin{bmatrix}\epsilon(1,1) \\ \epsilon(1,2) \\ \epsilon(1,3) \\ \vdots \\ \epsilon(20,20)  \end{bmatrix}$$
So $H$ is this matrix in the above equation, and it contains the PSF as a column.
Note that the $PSF_{\alpha,\beta}(r-r_0,c-c_0)$ and the image $d(r,c)$, make sense as two-dimensional objects, But they give you 400=20x20 separate equations, of the form in equation (2). You are just putting all those equations in 400 rows, and the columns in those equations make your matrix equation.
3

I have taken a y matrix (2D array)  ...  I don't know if it is correct to do this.

In my interpretation y should be a vector. Like I wrote explicitly in the above equation. The vector will have 20 x 20 entries.
Note that the vector/matrix notation only makes sense when you are applying linear algebra. E.g. when you solve the OLS problem. When you estimate $\alpha$ and $\beta$, the non-linear problem, then it is of little use.
4

Difference between deducing original signal and deducing parameters of PSF (Point Spread Function)

In the context of this problem you are estimating the parameters $a$, $b$, $c_0$, $r_0$, $\alpha$, and $\beta$. The first four parameters describing the original signal, and the latter two are parameters describing the PSF.
The difference is that deducing the original signal may turn the non-linear regression problem (when one of $c_0$, $r_0$, $\alpha$, and $\beta$ are unknown) into a linear regression problem (when $c_0$, $r_0$, $\alpha$, and $\beta$ are known).
