in question is

C:=2*sum(1/i^4, i=1..infinity);

OP's sum:
evalf[30](sum(sum(1/(i^4+j^4),i=1..infinity),j=1..infinity) + C);
      2.92637815841768603016360338185+0.*I

@Kitonum 

OK, I just wanted to have N as parameter, without thinking at speed.

It should have been:

eval(  'add'('add'('evalf'(a), i=1..N), j=1..N),  N=1000);

I understand op, nops as reflecting the logical representation of the expression.

@John Fredsted

Your question reduces to explain:

subs(a*b=c, 7*a*b); #  7*c
subs(a*b=c, I*a*b); #  I*a*b

The answer: distinct internal representation, as shown by

dismantle(7*a*b); dismantle(I*a*b);

namely, 7*a*b has two first level operands: 7 and a*b
but I*a*b has three (I, a, b).

The unpleasant fact is that op and nops do not reflect in this case the internal representation.

My impression is that the only problem is the limited capability of assume, as the next very simple examples show.

assume(p::integer, p>1, p<4);
assume(x<3,x>2);
about(p);
   Originally p, renamed p~:
     is assumed to be: AndProp(integer,RealRange(2,3))

[floor(x)-2, floor(p/2)-1, sin(floor(p/2)*Pi), cos(floor(p/2)*Pi)+1, sin(p-2)*(p-3)];



(should be all 0).

@Preben Alsholm 

I don't think that equality is the real issue here because <= is for sure a valid condition and p<=0,p>=0 gives p=0 (and about(p) confirms this).

@Preben Alsholm 

It seems that the assume database is not always inspected (correctly?):

[p, sin(p), sin(Pi*p)] assuming p=0;
  
    [p, sin(p),0]

@Markiyan Hirnyk 

int(x^(-ln(x)), x = 0 .. b);

also works (without assumptions)!

@JacquesC 

But most of these examples can be obtained in Maple 2015 as "one-liners", e.g.

Change(int(x^(-ln(x)),x=0..b), x=ln(t)) assuming b > 0;

Edit: I wanted x=exp(t), but even this strange change of variable works!

@Kitonum 

They are so easy to compose ... E.g.

J := Int((2+ln(x))/(sqrt((ln(x)+1)*x)*sqrt(1-x-x*ln(x))), x);

ListTools:-Search is not faster; it actually calls member. Just type:

showstat(ListTools:-Search);

The error can be easily solved for Maple 2015 (probably in 2016 it's ok):

with(IntegrationTools):
f:=t^2*exp(-t)/ln(t):
J:=Int(f,t=0..infinity):
Change(J,t=exp(s)):
Split(%,[-eps,eps]):
J1,J2:=op(1,%),op(3,%):
J1a:=subs(t=s, Change(J1,s=-t,t)) assuming eps>0:
Jeps:=normal(combine(J1a+J2)):
J:=eval(Jeps,eps=0);

evalf(J);
    2.153547692

But the following theoretical problem appears.

If f is a function defined in an interval (a,b) having a singularity  s in (a,b),
h is a C^2 diffeomorphism from (a,b) to (c,d) and h is used as a change of variables
int(f, a..b) ---->  int(F, c..d)

then the following situations are possible:

1. f is integrable in the sense of Cauchy principal value but F is not.
2. F is integrable in the sense of Cauchy principal value but f is not.

@Axel Vogt 

Just consider

Int(a/x + 1, x=0..1)

for a-->0.

Edit. The changes of variables are in the intervals [0,1-a] (where the integral is not improper)
and [1+a, oo) where the integrability is obvious. (0<a<1).

I also vote up yours because it's general (it appears in Fihtengolts, you probably know this excellent textbook).

@Markiyan Hirnyk

@Markiyan Hirnyk 

1. I know those theorems (I teach measure theory BTW).

2. I asked Maple to compute the limit of the integral (not the integral of the limit). So, in principle, you should ask Maple to base its computations.
[I admit that Maple probably did not check the conditions].

3. It's not necessary to have uniform convergence. It's not difficult to use the Lebesgue dominated convergence for the integrand of J1+J2 (which, BTW, is not your expression).