Replacing int with an approximate integrator and using LSSolve I find for the updated system this answer from LSSolve:
sol := [5.53989661895405, [d1 = -0.000305007492156158, d2 = 12.0947160055038, mu = -1894.38456773538]].
sol[1] is the sum of squares divided by 2. Thus Eq1^2+Eq2^2+Eq3^2 > 11 for the values found for d1, d2, and mu.
The signs of d1 and d2 are irrelevant since only their squares appear in the equations.
This doesn't prove that there is no solution, but it raises that suspicion.

An example of what I'm talking about is the following which uses a different integrator than I used for the result above:
 

RES_AP:=eval([Eq1, Eq2, Eq3],int=rcurry(Student:-Calculus1:-ApproximateInt,method=simpson,partition=51)):
sol:=Optimization:-LSSolve(RES_AP, variables={d1, d2, mu},iterationlimit=2000);
eval([Eq1, Eq2, Eq3],sol[2]);
IntegrationTools:-Expand(%);

The result here is
sol := [5.53989661896720, [d1 = -0.000973253903034116, d2 = -12.0947149137687, mu = -1894.38455370635]]
and the 3 residuals (i.e. Eq1,Eq2,Eq3)  are
[2.3533083414644, -2.3540885251539, -0.0005510445]

You could try a Riemann sum, i.e. replace the integral with one of the two sums

add( a2F[i]/w[i]*(w[i+1]-w[i]), i=1..N-1);
add( a2F[i]/w[i]*(w[i]-w[i-1]), i=2..N);

where N = 500.  As the final result you may choose the average of the two.
If the w-values are too far apart for that approch to make sense you could find a spline approximation for a2F using CurveFitting:-Spline. Then just integrate that approximation.

A home made example, where I construct my own data: W and F.
F is just sin applied to the values in W.

restart;
r:=rand(-0.01..0.01);
['r()'$5];
N1:=500:
W:=[seq(1..N1/10,0.1)] + ['r()'$N1-9]: 
N:=nops(W);
F:=sin~(W):
S1:=add(F[i]/W[i]*(W[i]-W[i-1]),i=2..N);
S2:=add(F[i]/W[i]*(W[i+1]-W[i]),i=1..N-1);
f:=CurveFitting:-Spline(W,F,w):
int(f/w,w=min(W)..max(W),numeric); # Takes a little while.
int(sin(w)/w,w=min(W)..max(W)); # For comparison.

In your worksheet Q is being defined as a function.
Try after the lines you have to write Q(t);
Then your integral is evaluated.
In Maple input (1D input and preferred by many) you could do:
 

restart;
P := t->r*Q[inf]/(2 + 2*cosh(-r*t + b));
int(P(s), s = 0 .. t);

Then you can try various simplifications.

You can do like this:
 

restart;

N:=1:M:=sqrt(N*(N+1)):N1:=1+N:w:=10:f:=1:                      
ini1:= x(0)=0.5,y(0)=0.5,z(0)=0;
dsys:={diff(z(t),t)=-(N1+M*cos(2*w*t))*z(t)-1+f*(x(t)+y(t)), diff(x(t),t)=-(N1-I*w-2*M*exp(-2*I*w*t))*x(t)-f*(N1+(z(t)))-2*f*M*exp(2*I*w*t),diff(y(t),t)=-(N1+I*w-2*M*exp(2*I*w*t))*y(t)-f*(N1+(z(t)))-2*f*M*exp(-2*I*w*t)};
### The crucial difference:
zd:=subs(dsys,diff(z(t),t));
###                          
res:=dsolve(dsys union {x(0)=0.5,y(0)=0.5,z(0)=0},numeric,output=listprocedure);

tit:=sprintf("F=%g,N=%g",f,N);

plots[odeplot](res,[[t,zd]],0..6,axes=boxed,color=blue,tickmarks=[3, 4],thickness=2,title=tit,size=[700,default],numpoints=1000);

##TRY:
combinat:-choose([$(0..n-1)], k);
add(mul~(combinat:-choose([$(0..n-1)], k))) ;

You aren't going to get explicit solutions for this equation. But you do get some good information from an initial symbolic treatment.
After that use numerical solution.
Notice that because a'(t) is squared dsolve will first have to find the 2 solutions for a'(t) and then solve those, in this case using separation of variables. The first result that pops up besides the 2 mentioned corresponds to the exceptional case when separating variables: that a denominator can be zero.
 

restart;
ode:=3*a(t)-a(t)^2-3*z*diff(a(t), t)^2+c*a(t)+w/a(t) = 0; # c*k replaced by just c
dsolve(ode); # 3 results
value([%]); #Real large so forget it.
sol:=dsolve(ode,implicit); # Neater
## Solving for the derivative:
odes:=solve(ode,{diff(a(t),t)});
ode1,ode2:=op(odes[1]),op(odes[2]);
dsolve(ode1); # implicitly given
## Finding equilibria:
RHS:=subs(a(t)=a, rhs(ode1));
E:=solve(RHS=0,a); #Equilibria: Same as solving sol[1] for a(t)
nops([E]); # 3 as expected, but not all are real
evalf(eval(E,{c=1,w=2,z=3})); # Example
e1:=%[1]; # The real one
eval(sol[1],{a(t)=e1,c=1,w=2,z=3}); # basically 0=0 as it should be
## 
## Numerical solution using the parameters option:
res:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z]);
res(parameters=[1,2,3]); # Setting parameters of choice.
plots:-odeplot(res,[t,a(t)],-1..10,caption=typeset("Reaching equilibrium ",e1," in finite time"));
## The equilibrium is reached in finite time because of a sqrt(e1-a(t)) in ode1.
## We could find the time of arrival to e1 this way, if we want:
resE:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z],events=[[sol[1],halt]]);
resE(parameters=[1,2,3]);
resE(8); # Time 5.9388290

Just name the inputs, as in
 

ode := 0.205*diff(u(t), t, t) + 0.53776*diff(u(t), t) + 7.49*u(t);
ics := u(0) = 0.6039948504, D(u)(0) = 0.2696963684;
dsolve([ode, ics]);

Begin by trying Y in this modified version where I only give the 2 procedures. Everything else is the same:
 

# Procedure to return values for specific time and parameter values.
u := proc( t, a, b ) 
	sol( 'parameters' = [ ':-a' = a, ':-b' = b ] ):
	try
		return eval( y( ':-t' ), sol( t ) ):  ### Changed
	catch:
		return 10000:
    end try:
end proc:


# Procedure to return objective function.
r := proc( a, b )
	local A, t, p, q, x, y:
	A := [ seq( u( t, a, b ), t=T ) ] -~ Y:  ### Changed
	return foldl( (p,q) -> p + q^2, 0, op(A) ):
end proc:

With X instead of Y Minimize cannot find an improved version.
A modified version allowing for several lists X,Y,Z, .. would consist in just replacing the zip part as follows:
 

`[]`~(X,Y) ## replacement for zip( (x,y) -> [x,y], X, Y )

This would work with no other changes in the original code for the 2 procedures. In the X,Y case the full line would be:
 

A := ListTools:-Flatten( [ seq( u( t, a, b ), t=T ) ] -~  `[]`~(X,Y) );

If only using one of them just remove the other, as in  `[]`~(Y).
One minor thing:
The return value in the catch error case ought to be of the same type as in the no error type:
i.e. [10000,10000] when using X,Y and [10000] when only using Y.

A note about the failure of Minimize when using X.
In the exact result for the system we see that x(t) only depends on the product a*b, whereas y(t) depends on both a and b independently.
This is going to be a problem for any numeric solver:
Try this version tailored to that situation:
 

ux := proc( t, ab ) #name ab. You could use any name of course
	sol( 'parameters' = [ ':-a' = ab, ':-b' = 1 ] ): # Arbitrarily setting b = 1.
	try
		return eval( [ x( ':-t' ) ], sol( t ) ):  
	catch:
		return [10000]:
    end try:
end proc:


# Procedure to return objective function.
rx := proc( ab ) #one input ab
	local A, t, p, q, x, y:
	A := ListTools:-Flatten( [ seq( ux( t, ab ), t=T ) ] -~  `[]`~(X) );  
	return foldl( (p,q) -> p + q^2, 0, op(A) ):
end proc:

rx(3);
# Find parameter values for bet fit.
Optimization:-Minimize( 'rx'(ab), ab=-1..1); # One variable ab

If you eliminate i(t) as done below your problem only depends on 3 parameters instead of 5.
That ought to help no matter what the method is.
 

restart;

Digits:=15:
eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t);
eq2:=L*diff(i(t),t)+R*i(t)=V-K*diff(theta(t),t);

ICs:= theta(0)=0, D(theta)(0) = 0, i(0) = 0;
####
M:=ImportMatrix("F:/MapleDiverse/servo_open_ident_1_A.csv",source=csv):
M[1..10,..]; #Jump in V (col. 2) at start
M[-10..-1,..];
plot(M[..,1..2],labels=[t,V]); 
#### For V I shall use this:
V1:=Statistics:-NonlinearFit(a+b*sin(c*t+d),M[..,1..2],t,initialvalues={a=0,b=0.5,c=1/50,d=0});
plot(V1,t=0..50); p1:=%:
plot(M[..,1..2]); p2:=%:
plots:-display(p1,p2);
w1:=unapply(V1,t)~(M[..,1]):
w:=M[..,2]:
dw:=w-w1:
eval(V1,t=0);
dw[1];
(max,min)(dw[2..]); # Not bad
####
plot(M[..,[1,3]],size=[1000,default]); pd:=%:  # For later use kept in pd.
plot(M[..,[1,4]],size=[1000,default]); # Isn't used in the following
params:={J= 6.9347e-005, b = 2.2480e-004, K= 0.0094, R=0.2124, L=0.0015}; # Supplied by OP
####From eq1 and ICs follow that (D@@2)(theta)(0)=0. Will be used below.
####Eliminating i(t) by using integrating factor exp(R/L*t) used to solve first order linear odes:
S:=Diff(exp(R/L*t)*i(t),t)=exp(R/L*t)*rhs(eq2)/L;
##Operating on eq1:
map(Diff,exp(R/L*t)*convert(eq1,Diff)/K,t);
##Elimination of i(t):
subs(S,%);
ode:=value(%);
ode1:=L*J*op(solve(ode,{diff(theta(t),t,t,t)}));
ode2:=collect(ode1,diff,factor);
##We see that this third order ode depends on the following groups only: 
TR:={c1=(J*R+L*b)/(L*J), c2=(K^2+R*b)/(L*J), c3= K/(L*J)};
cparams:=eval(TR,params); #Op's params translated to c1,c2,c3
sbs:=solve(TR,{J,K,b});
ode3:=subs(sbs,ode2);
solve(ode3,{diff(theta(t),t,t,t)});
ODE:=eval(op(%),V=V1); # The 3rd order ode for theta(t).
##The initial conditions:
ics:=theta(0) = 0,D(theta)(0) = 0,(D@@2)(theta)(0)=0;
res:=dsolve({ODE,ics},numeric,parameters=[c1,c2,c3],output=listprocedure,maxfun=10^6);
th:=subs(res,theta(t));
## Q sets the parameters and computes the sum of squares of the residuals:
Q:=proc(c1,c2,c3) option remember; local j,ssq;
  if not [c1,c2,c3]::list(realcons) then return 'procname(_passed)' end if;
  try
    res(parameters=[c1,c2,c3]);
    ssq:=add( (th(M[j,1])-M[j,3])^2,j=1..5001)
  catch:
    ssq:=10^15
  end try;
  ssq
end proc:
##
t0:=time[real]():
RES:=Optimization:-NLPSolve(Q(c1,c2,c3),initialpoint=cparams);
time[real]()-t0; # 640s in my case.
##Finished with the error: No improved point could be found.
##But let us look at the remember table for Q:
T:=op(4,eval(Q)):
LL:=[entries(T,pairs)]:
LL:=remove(hastype,LL,{function,HFloat(undefined)}):
LS:=sort(LL,key=rhs):
LS[1];
LS[-1]; # A negative value for c1
res(parameters=([c1,c2,c3]=~[lhs(LS[1])]));
plots:-odeplot(res,[t,theta(t)],0..50,color=blue); pm:=%:
plots:-display(pd,pm,legend=[Data,Model]); #Not bad

The final plot of theta:

The c-values found were: [c1 = 0.07085248546343231, c2 = 1072.44657764141, c3 = 91326.7791941768]
(the result of [c1,c2,c3]=~[lhs(LS[1])]; )

Why not begin by solving this simple linear system exactly?
 

restart;
J:= 0.01;b:=0.1;K:= 0.01;R:=1;L:=1;
eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t);
eq2:=L*diff(i(t),t)+R*i(t)=V-K*diff(theta(t),t);

ICs:= theta(0)=0, D(theta)(0) = 0, D(theta)(0) = 0, i(0) = 0:

res:=dsolve({eq1,eq2,ICs});

Then take it from there?

I suppose you might want to remove all equations of the form a = a ?
Then this would do it:
 

remove(x->(lhs-rhs)(x)=0,EQI__4);

or shorter:
 

remove( lhs-rhs = 0, EQI__4);

The result of splitting and a change of variable below seems alright, but the first result without splitting is I'm afraid a result of carelessness.
 

restart;
J := Int(cos(2*x)/(1+2*sin(3*x)^2), x = 0 .. Pi);
J1:=combine(J);
value(J1); # returns unevaluated (except with int instead of Int)
IntegrationTools:-Change(J1,cos(2*x)=z); # result right, but is it due to carelessness?
S:=IntegrationTools:-Split(J1,[Pi/4,Pi/2,3*Pi/4]);
L:=convert(S,list); # To see the steps
map(IntegrationTools:-Change,L,cos(2*x)=z);
convert(%,`+`);

Your system satisfies the requirements of uniqueness for an initial value problem, and that is what you got.
So since r[1] = V = q = 0 is a solution, it is the unique solution.

@siddikmaple Try Optimization:-LSSolve instead:
 

CSTR:={p[b], p[s], p[0, b], p[0, c], p[0, p], p[0, s], p[1, b], p[1, c], p[1, p], p[1, s], p[2, b], p[2, c], p[2, p], p[2, s], p[3, b], p[3, c], p[3, p], p[3, s], r[0], r[1], r[2], r[3], t[0], t[1], t[2], t[3]}=~(0..1);
##
residuals:=(lhs-rhs)~([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14,e15,e16,e17,e18,e19,e20, e21,e22, e23, e24, e25, e26]);
##
ansLS:=Optimization:-LSSolve( residuals,op(CSTR),iterationlimit=10000);

The answer given is:

Continue with using that result as an initialpoint in SolveEquations:
 

DirectSearch:-SolveEquations(residuals,CSTR,initialpoint=ansLS[2]);

The answer is the same. Notice that LSSolve returns half of the sum of squares and gets:
0.563125933394723333e-1
whereas SolveEquations returns the sum of squares and finds:
.112625186678945
So there is basic agreement.

@siddikmaple You cannot use [ ] or { } instead of parentheses ( ) .
[ ] is used for lists and also for selection. { } is used for sets.
You have 4 cases , where [ ] is used instead of ( ), and 5 cases where { } is used instead of ( ):

sys:={e1, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e2, e20, e21, e22, e23, e24, e25, e26, e3, e4, e5, e6, e7, e8, e9}:
nops(indets(sys,list)); # 4
nops(indets(sys,set)); # 5

This must be corrected before having any chance of getting a result out of fsolve.