How is it possible to collect NLexpression in a loop?

Would you please help me that how I can collect some NLexpressions in a loop? I Want to keep k8 and k9 for all i=1:10 as scenarios. It means in the end of the loop, Q equal to collection of k8 and k9 under each scenario ( i ). I couldn't define a matrix and put each pair of k8 and k9 in that as an element. with considering Q the code doesn't work as well. Many thanks for your kindly help.

using JuMP,CPUTime, Distributions, Ipopt,Juniper,Cplex
n1=1; #the least of scenarios
N=4; #number of scenarios
M=20; #number of sampling
landa=0.01;
E=0.05
T0=0;
T1=2;
T2=2;
gam2=1; gam1=1;
a1=0.5; a2=0.1; a3=50; ap=25;
c0=10;
Zn=zeros(N, 4)
Q=0;
for i in n1:N
    C1=rand(100:100:300);
    sig=rand(0.5:0.5:2);
    f(x) = cdf(Normal(0, 1), x);
    #---------------------------------------------------------------------------
    ALT= Model(optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0),

           "mip_solver"=>optimizer_with_attributes(Cplex.Optimizer, "logLevel" => 0),"registered_functions" =>[Juniper.register( :f, 1, f; autodiff = true)])

           );

    # variables-----------------------------------------------------------------
    JuMP.register(ALT, :f, 1, f; autodiff = true);
    @variable(ALT, h >= 0.001);
    @variable(ALT, L >= 0.000001);
    @variable(ALT, n>=2, Int);

    #---------------------------------------------------------------------------

    @NLexpression(ALT,k1,h/(1-f(L-sig*sqrt(n))+f(-L - sig*sqrt(n)))); # HARL1

    @NLexpression(ALT,k2,(1-(1+landa*h)*exp(-landa*h))/(landa*(1-exp(-landa*h)))); #to

    @NLexpression(ALT,k3,E*n+T1*gam1+T2*gam2);
    
    @NLexpression(ALT,k8,(C1*(k1-k2+k3)));# depend on scenario
    
    @NLexpression(ALT,k9,(((a1+a2*n)/h)*(k1)));#depend on scenario

    Q=Q+k8+k9;
    #-----------------------------------------------------------------------
end

You have a couple of issues in your code.

• Why are you creating a new model every loop? You cannot aggregate expressions across models like you do with Q=Q+k8+k9
• You cannot add nonlinear expressions like Q=Q+k8+k9. All nonlinear expressions must occur inside macros

In general, you are not limited to the specific syntax of JuMP. You can use any Julia data structures to help. In the code below, I just push expressions into a vec

I haven't tested, so there might be typos etc, but this should point you in the right direction:

using JuMP, Distributions
E, landa, T1, T2, gam1, gam2, a1, a2 = 0.05, 0.01, 2, 2, 1, 1, 0.5, 0.1
ALT = Model()
f(x) = cdf(Normal(0, 1), x)
JuMP.register(ALT, :f, 1, f; autodiff = true)
@variable(ALT, h >= 0.001)
@variable(ALT, L >= 0.000001)
@variable(ALT, n >= 2, Int)
k8, k9 = Any[], Any[]
for i in 1:4
    C1 = rand(100:100:300)
    sig = rand(0.5:0.5:2)
    k1 = @NLexpression(ALT, h / (1 - f(L - sig * sqrt(n)) + f(-L - sig * sqrt(n))))
    k2 = @NLexpression(ALT, (1 - (1 + landa * h) * exp(-landa * h)) / (landa * (1 - exp(-landa * h))))
    k3 = @NLexpression(ALT, E * n + T1 * gam1 + T2 * gam2)
    push!(k8, @NLexpression(ALT, C1 * (k1 - k2 + k3))
    push!(k9, @NLexpression(ALT, k9, (a1 + a2 * n) / h * k1)
end
Q = @NLexpression(ALT, sum(k for k in k8)  + sum(k for k in k9))