Hello Julia community. I am working on a heat equation of 1D heat rod on which I want to apply the Model Predictive Controller. I have pasted my edited code with this message. I wanted some help to how can I solve the errors. Some of the suggestions that I have obtained are as follows:

(1) I need to change the position as to which is the best place to define the model.

(2) Which coefficients do my MPC algorithm need and where do I implement them? And what are suitable values for them?

e.g. Sampling time, MPC horizon, final simulation time, etc…

(3) ow do I implement the system dynamics / Euler step z[k+1] = z[k] + ΔT * f(z[k], u[k]) correctly?

Does not work: mech_osc(θ,u) = α/(Δx^2)* θ + 2/(c * ρ * Δx) * u

How should I change my 2 “for” loops such that the part “@constraint” gets removed.

Any small help will be great for me. Looking forward to your suggestions.

Thank you!!

```
using JuMP, Ipopt, LinearAlgebra
mod = Model(Ipopt.Optimizer)
set_optimizer_attribute(mod, "max_iter", 100)
set_optimizer_attribute(mod, "mumps_mem_percent", 500)
#A JuMP model is initialized
#The parameters of the 1-D rod
L = 0.1 #Length of rod
λ = 45.0 # Conductivity
c = 460.0 # Specific heat capacitivity
ρ = 7800.0 # Material density
α = λ / (c*ρ) # Diffusivity
#Number of spatial discretization points and length of each
N = 11 #Discretization points / nodes
Δx = L / (N-1) # Δx = x[i+1] - x[i]
#Integration from t0 = 0 to tf
tf = 900.0 #Final time: 900 seconds = 15 minutes
Δt = 1.0 #Sampling period
k = round(Int,tf / Δt); #k = number of steps
#Initial and reference temperatures in Kelvin
θinit = 273.0
θref = 500.0
#bounded input values[u_min,u_max]
u_min = 0.0
u_max = 20e5;
function mpc_emo(θinit)
#Decision variables of the optimization problem
@variables mod begin
θ[1:N,1:k]; # temperature at x[1:N] at time step k
u_min <= u[1:N, 1:k - 1] <= u_max; # bounded constrained Input u
y[1:k]; # Output y
end
#Initial values
@constraint(mod, θ[:,1] .== θinit)
@constraint(mod, y[1] .== θinit)
#System dynamics
for j in 1:(k-1)
for i in 2:N-1
@constraint(mod,θ[i,j + 1] == θ[i,j] + Δt * α/( Δx^2)(θ[i-1,j] - 2θ[i,j] + θ[i+1,j]) ) # Diffusion at inner grid points
end
@constraint(mod, θ[1,j + 1] == (θ[1,j] + Δt*( α/(Δx^2) * (-2*θ[1,j] + 2*θ[2,j]) + 2/(c * ρ * Δx)*u[j]))) # Left side: heat input
@constraint(mod, θ[end,j + 1] == ( θ[end,j] + Δt * α/(Δx^2) * (2*θ[end-1,j] - 2*θ[end,j]) )) # Right side: isolated side
@constraint(mod,y[j + 1] == θ[end,j + 1]) # Measurement at left side
end
q_w = 1000; # Weighing coefficient for errors (or states)
r_w = 1e-4; # Weighing coefficient for input signals
# e(t) = r(t) - y(t)
# err = q * sum( e(t_i)^2 ) with index i ∈ (1,k-1)
err = @NLexpression(mod, sum(q_w * (θref - y[j])^2 for j in 1:k-1) )
#in_err = r * sum( u(t_i)^2)
in_err = @NLexpression(mod, sum(r_w * u[j] ^ 2 for j in 1:k-1) )
# Cost function
J = @NLexpression(mod,0.5 * Δt * (err + in_err))
#defining the cost function
@NLobjective(mod, Min, J)
#optimizes the model
optimize!(mod)
return JuMP.value.(θ), JuMP.value.(u)
end
mech_osc(θ,u) = α/(Δx^2)* θ + 2/(c * ρ * Δx) * u
θpos = rand(11) # Recent states
θ_hist = zeros(11,0) # History of states
for i in 1:50
# Run the MPC control optimization
θ_plan, u_plan = mpc_emo(θpos)
# Save states
θ_hist = hcat(θ_hist, θ_plan[:,1])
θpos = θpos + Δt * mech_osc(θpos,u_plan[:,1])
end
using Plots
plot(θ_hist')
```