This has nothing to do with you loading packages Statistics and SignalProcessing (even if this can generate some confusion)
The fact that Statistics:-AutoCorrelation gives a wrong result, is more profound than that.

The formula used in this function is basically wrong and based on a kind of generalisation of the statistical correlation between two vectors.
The best proof this formula is wrong, IMO, comes from the computation of the Statistics:-CrossCorrelation : CrossCor(V, V) = AutoCor(V) up to a shifting of N/2 along the horizontal axis.

You will see in the attached file that (up to this x-shifting) CrosCor(V, V) is correct as AutoCor(V) is not.
Comparing the formulas used, CrosCor(V, V) operates on the raw vector V, while AutoCor(V) operates on the centered W=V-Mean(V) (mu= Mean(V) hereafter)
But Autocor(W)[k] = CrossCor((V-mu), (V-mu))[k] ~ CrossCor(V, V)[k]- mu^2*(N-k) (which becomes negative for k large enough)

Here are a lot of comparisons of different wys to compute the autocorrelation function.

autoC_mmcdara.mw

This way

Download plotsetup_ps.mw

I can't visualize ps file on my Imac and online ps document viewer don't seem compatible with the ps file generated by Maple.

So I saved with jpeg format (plotsetup(jpeg, ...)) used an online converter (https://anyconv.com/fr/convertisseur-de-jpg-en-ps/) to generate the desired ps file.
The result can't be loaded on this site bus is in all points identical to the image below.

Here is the ps file (I added a pdf extension in order to upload it).
AnyConv.com_plot.ps.pdf
Remove the pdf extension and visualize the ps file.

EDO b5 and b6 both contain the unknown function f (= f(x)), which is not allowed if you use the option numeric.

Remedy: instantiate f and apply dsolve; for instance

c5 := dsolve({bcs1, eval(b5, f(x) = x)}, numeric); c5(0); 
c6 := dsolve({bcs1, eval(b6, f(x) = x)}, numeric); c6(0);

Same comments as  @dharr, one fixed the "parameters" problem it remains the one of Q(t).
Covid19_Simulation_mmcdara.mw

What happens if you use the first strategy described at the end of this worksheet:

# First strategy
# define Q(t) (notional example)
MyQ := Q(t) = t:

# remove tee last IC
init_conds := subsop(-1=NULL, [init_conds]):

# form the new system
sys    := {Asymp, Exp, Immun, Inf, Quar, Suscep, Vacc, init_conds[]}:
sys    := eval(sys,MyQ):
params := convert(indets(sys, name) minus {t}, list):

# dsolve it
sol    := dsolve(sys, numeric, parameters = params, method = rkf45):

Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

The reason is likely that you have two different equations of the same form:
                            diff(I(t), t) = xxx  and  diff(I(t), t) = yyy  with xxx <> yyy
So you need to remove one or write something like

diff(I(t), t) = piecewise(condition, xxx, yyy)

For instance, removing the equation Inf enables dsolve to find a solution:

sys    := {Asymp, Exp, Immun, Quar, Suscep, Vacc, init_conds[]}:
sys    := eval(sys,MyQ):
params := convert(indets(sys, name) minus {t}, list):
sol    := dsolve(sys, numeric, parameters = params, method = rkf45)
                  proc(x_rkf45)  ...  end;

As I don't want to interfere with @Carl Love I'll just give you a quick reply.
 

P := (p, q) -> X -> op(0, X)(subsop(p=[op(X)][q], q=[op(X)][p], [op(X)])[]):
P(2, 3)(X(i, j, k));
                           X(i, k, j)
P(3, 2)(%);
                           X(i, j, k)
P(2, 3)(U("a", 3, z));
                          U("a", z, 3)
P(3, 2)(%);
                          U("a", 3, z)

A sketch of a more general permutation than a 2-by-2 exchange

P := (Before, After) -> X -> op(0, X)(subsop((Before =~ [op(X)][After])[], [op(X)])[]):
P([1, 2, 4], [4, 1, 2])(F(p, q, r, s, t))

                        F(s, p, r, q, t)

Here the detail of the exchanges is

p, q, r, s, t;  # original sequence
s, q, r, s, t;  # op 1 receives the value of op 4 from the original sequence
s, p, r, s, t;  # op 2 receives the value of op 1 from the original sequence
s, p, r, q, t;  # op 4 receives the value of op 2 from the original sequence

Thus subsop does the replacements all at once by refering to the original sequence.
Maybe you could want these replacements made iteratively, in which case you can get 

p, q, r, s, t; # original sequence
s, q, r, s, t; # op 1 receives the value of op 4 from the original sequence
s, s, r, s, t; # op 2 receives the value of op 1 from the previous sequence
s, s, r, s, t; # op 4 receives the value of op 2 from the previous sequence

of what @dharr already did.

The first one seems closer to to what you have been asked.
The second one is more concise as it uses GraphTheory features, but is less close to your step by step request.

Download Graph_mmcdara.mw

No procedure is given so that you can see the results of the different steps.

Here are two codes, one for multivariate PLS, the other for PCR.

As I told you from the office (sand15), I extracted these codes from a larger application which contains it's own (proprietary) documentation.
So it will be probably quite difficut to understand what the two codes return (in particular for PLS) if you are not already familiar with the two factor analysis they implement.
I've tried to give a minimum of informations and references, specifically for PLS. Unfortunately (maybe?) as I'm French these references are a course and a text book both written in french.

I don't know what is your knowledge in factor analysis, but all these methods, and even the "simplest" PCA (see MAPLE help pages), require a certain amount of practice for their outputs to be correctly interpreted.
This is all the more tru for PLS.
If you have a correct understanding of regression and PCA, you should not have any difficulty to understand PCR outputs (which I named ACPVI in my worksheet) for it is nothing but a regression on the synthetic variables a PCA produces.

This first file contains PLS and PCR algorithms plus a test case named "Linnerud"
see here for a short description of these data  https://www.rdocumentation.org/packages/ropls/versions/1.4.2/topics/linnerud
PLS+PCR_Linnerund.mw

This second file contains PLS algorithm plus a test case named "Russett"
see here for a short description of these data  https://search.r-project.org/CRAN/refmans/RGCCA/html/Russett.html
PLS_Russett.mw

A in all factor analysismethods, more important than the algorith itself are the many plots or matrices you have to display to unfold the analysis correctly.
The variety of these "analysis tools" is that large that only a few are provided in these two worksheets.

So, please consider the PLS I delivered you as a "MINIMAL" version of the method. Writting a complete  and operational implementation of PLS in Maple is equivalent to develop a complete dedicated package and would take a week or so, not speeaking of the help pages.
Go get an idea of that look here the content of the R pls package:
https://cran.r-project.org/web/packages/pls/pls.pdf

If you just want to have indentation (for instance for "for...do " or a "if..then..else" structure) to verfy your code, you can use NOTEPAD++ :

• export your worksheet as mpl format.
• open the mpll within NOTEPAD++
  (ice on the cake: the Maple code will be displayed using syntactic coloring)
  here is an old thread https://www.mapleprimes.com/posts/129552-Maple-In-All-Colours-Of-Rainbow-
   

If you want to have automatic indentation like "in other languages" while writting your code in the worksheet, my answer is "I don't know if it'is possible".
Wait for someone else's answer

Firstly remove these unavoided "*" after display, draw and textplot.

Secondly, it's not 

draw[A(color = black, symbol = solidcircle, symbolsize = 6), 
B(color = black, symbol = solidcircle, symbolsize = 6), 
C(color = black, symbol = solidcircle, symbolsize = 6), 
ABC(color = blue)]

but

draw([A(color = black, symbol = solidcircle, symbolsize = 6), 
B(color = black, symbol = solidcircle, symbolsize = 6), 
C(color = black, symbol = solidcircle, symbolsize = 6), 
ABC(color = blue)])

use the big green up-arrow to load your worksheet.

Nevertheless the message is rather clear: you use a variable named "v" that Maple doesn't know what to do with....

The reason is also rather clear: you claimed you need "help in solving the following system of partial differential equations" but you hvaveonly differential equations.
Alllf the unknown equations a(t, v), c(t, v) ... appear in expressions like diff(a(t, v), t), diff(c(t, v), t) ...
Then "v" is something whish is not (from the point of view of dsolve(..., numeric) ) a dependent variable, nor an independent one.

Lookk at this piece of code

Download WhereYourErrorComesFrom.mw

I suggest you correct yourself your own worksheet and come back to this same thread if you still have difficulties (which is very likely given the equations eq1..eq9 you wrote).

Several points more or less already emphasized by @dharr and @Carl Love

• ctrs8 lacks a parenthesis somewhere at the end:

  # is it 
  0.1749271137e-1*(....)^2)-1)
  # or
  0.1749271137e-1*(....))^2))-1 
  # ?

• ATS contains 10 indeterminates and you want to Minimize/Maximise a function of NN and MM

  indets(ATS(NN, MM), name);
  numelems(%);
     {MM, NN, H[1], H[2], H[4], L[1], L[4], L[7], P[X1], P[XF]}
                                 10
  

  You need to give H[1], H[2], H[4], L[1], L[4], L[7], P[X1], P[XF] numerical values before trying to run 

  NLPSolve(ATS, {cstr8}, bounds, initialpoint = [15.90, 18.90], assume = nonnegative)

• Your description of cstr8 is almost unbereable without loadind your wotksheet.
  Here is something clearer

  cstr8 := proc(NN, MM) 
  0.1749271137e-1*
  (
    588.000000000000+72.0000000000000
    *
    (-3.03333333+.333333333333333*MM)
    *
    (MM-17.80)
    +
    72.0000000000000
    *
    (.333333333333333*NN-3.5666666666667)
    *
    (17.00-1.*NN)
    +
    24.50000000*
    (2.90+4.666666667*tan(arctan(.166666666666667*MM-2.96666666666667)))
    *
    (18.00+5.627472876000
      *
      sin(-1.570796327+arctan(.166666666666667*MM-2.96666666666667))
      /
      sin(-1.745329252+arctan(.166666666666667*MM-2.96666666666667))
    )
    *
    cos(-.5235987756000+arctan(.166666666666667*MM-2.96666666666667))^2
    *
    sin(
      .3490658504000
      +
      arctan(
        .166666666666667*MM
        -2.9666756000
        +arctan(.166666666666667*MM-2.96666666666667)
      )^2
    )
    -
    1
  )
  end proc:
  

  (please validate the position of the last parenthesis)
  Using map(identify, ...)  here is an even more consise form of cstr8 (please validate the position ofthe last parenthesis)

  cstr8 := proc(NN, MM) 
  6/143 * (587+(72*(-91/30+(1/3)*MM))*(MM-89/5)+(72*((1/3)*NN-107/30))*(17-NN)+(49/2*(-197/18+(7/9)*MM))*(18+5.627472876*sin(-(1/2)*Pi+arctan((1/6)*MM-89/30))/sin(-(5/9)*Pi+arctan((1/6)*MM-89/30)))*cos(-(1/6)*Pi+arctan((1/6)*MM-89/30))^2*sin((1/9)*Pi+arctan((1/6)*MM-2.9666756+arctan((1/6)*MM-89/30))^2))
  end proc;
  

• You miss 4 extra constraints representing NN and MM bounds:

  cstr1 := proc(NN, MM) 14.90-NN end proc:
  cstr2 := proc(NN, MM) NN-16.90 end proc:
  cstr3 := proc(NN, MM) 17.90-MM end proc:
  cstr4 := proc(NN, MM) MM-19.90 end proc:
  

Given all this here is the correct syntax to get a solution (which depends ont the values you gave to H[1], H[2], H[4], L[1], L[4], L[7], P[X1], P[XF])

Optimization:-NLPSolve(
  fATS
  , {cstr1, cstr2, cstr3, cstr4} #, cstr8}
  , initialpoint=[15.9, 18.9]
)

Full code here (it's up to you to validate cstr8)
(only an excerpt is displayed to whow the instanciation of the the parameters and the solution)

Download NLPSolve_mmcdara.mw

I agree one hundred percent with @tomleslie "it is simpler to generate the data used by the plot".

However if, for some reason you truly want to "extract" the points (in fact segemnts of straight line) that define each contour displayed by contourplot, here is a way:

# assuming this is your plot
p := contourplot(x^2+y^2,x=-1..1, y=-1..1, contours=2):
display(p);

# then levelcurves contains a list of 2 lists (because I set contours=2) of segments
C := select(has, [op(p)], CURVES):
levelcurves := map(u -> select(type, [op(u)], list), C):

One problem can be that each list, even if the contour it represents is defined by a single piece (in my exemple the inner circle), 
Even if contour is connected (made of a single curve, such as the inner circle in my example) one problem can the segments that the corresponding list of levelcurves contains, are not contiguous in the following sense:

• Let L the list from levelcurves which describes the inner circle.
• Let N the number of segments L contains.
• Let S[n] [departure (1 <= n < N) one such segment. 
• Then arrival(S[n]) is not necessarily equal to departure(S[n+1])

This is detailed in the attached worksheet.
This worksheet provides a piece of code to rearrange a "levelcurve" in a list where arrival(S[n]) = departure(S[n+1]) for each n.

Contiguous_Curve_from_Contourplot.mw

When the "levelcurve" is not connected this kind of sorting becomes more abiguous (one must define and order on the segments) and more complex (the sorting must be made for each connected component).
Note that this same "unordered" representation of a curve also appears with plots:-intersectplot.

@dharr already gave you some informations but I understand that what you are interested in is to asess the variance matrix of the parameters?

This is possible with Statistics:-Fit for this matrix can be build from simple linear algebra relations.
For non linear fit there is no simple way to assess this matrix and, as it is always the cas in this kind of situations, the "solution" has name "resampling".

The basic idea of resampling is to simulate what we would have obtained had we collected data from other sets of experiments.
For instance R people realize independently the same set of N experiments and, still independently, asses the parameter vector.
This give a set of R estimations of the parameter vector we would obtain if we realize an infinite set of  experiments.
But the important point is that we can now build an estimation of the covariance matrix of the parameter vectorby using these R "replicates".

This is explained in detail in the file below for yourfirst model.
 

Download NonlinearFit_CorrelationMatrix.mw

A remark: In the future, could you load your worksheet using the green up arrow? 
This will simplify the work of those who would like to answer you.
Thanks in advance for all of us.

restart
with(plots):
Fig := proc(t) 
local a, b, P, Q, N, R, TG, x0, y0, p1, p2, p3, po, tp, sol; 
a := 11; 
b := 7; 
R := sqrt(a^2 + b^2); 
P := [R*sin(t), R*cos(t)];
x0 := P[1];
y0 := P[2]; 
TG := (a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*y0*x0*(x - x0)*(y - y0) = 0; 
p1 := implicitplot(x^2/a^2 + y^2/b^2 - 1, x = -11 .. 11, y = -7 .. 7, color = blue); 
p2 := implicitplot(x^2 + y^2 - a^2 - b^2, x = -15 .. 15, y = -15 .. 15, color = blue); 
p3 := implicitplot(TG, x = -15 .. 15, y = -15 .. 15, color = red); 

sol := solve(({x^2/a^2 + y^2/b^2 - 1, TG}), {x, y}, explicit); 
if sol = NULL then
  sol := solve(evalf({x^2/a^2 + y^2/b^2 - 1, TG}), {x, y}, explicit); 
  sol := remove(has, {sol}, I);
end if:

if sol <> {} then 
  Q := [subs(sol[1], x), subs(sol[1], y)]; 
  N := [subs(sol[2], x), subs(sol[2], y)]; 
  po := plot([P, Q, N], style = point, symbolsize = 15, symbol = solidcircle, color = red); 
  tp := textplot([[P[], "P"], [Q[], "Q"], [N[], "N"]], 'align' = {'above', 'left'}); display([p1, p2, p3, po, tp], scaling = constrained); 
  return display(po, tp); #  is it this that you want to plot ???
end if:
end proc:

nFig := 60:
display(seq(Fig(2*Pi*i/nFig), i=0..nFig), insequence = true):

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/A_solution.mw .
 

Download A_solution.mw

For Windows XT to Windows 10 and for Maple 2015 to Maple 2020 :

MyWS contains thename of the active WorkSheet IF and ONLY IF you have ONLY ONE MAPLE SESSION OPENED.

Be carefull, as I use a Windows version purchased from the French distributor, some informations are in French.
For instance, the word "window" translates to "fenêtre" in French and it is this word that appears in the return from the tasklist command. 
On the line tagged ####### you will have to change the field  "fenˆtre: " (note that Windows doesn't spell correctly the word "fenêtre") by the correct field.
T see what the correct field must be replace StringTools:-Squeeze(MyTask): by StringTools:-Squeeze(MyTask);
 

Download Get_MW_name.mw

For Mac OSX

I was never capable to get the name of the active WorkSheet on my Mac

For Linux

Never tried to do it.