src/IPM/ipmdata.jl
changed.
src/IPM/MPC/MPC.jl
changed.
src/IPM/HSD/step.jl
changed.
src/IPM/MPC/step.jl
changed.
src/IPM/HSD/HSD.jl
changed.
130 | 130 | c0 = pb.obj0 |
|
131 | 131 | if !pb.objsense |
|
132 | 132 | # Flip objective for maximization problem |
|
133 | - | c .= -c |
|
133 | + | c .= .-c |
|
134 | 134 | c0 = -c0 |
|
135 | 135 | end |
|
136 | 136 |
170 | 170 | u = [pb.uvar; uslack] |
|
171 | 171 | ||
172 | 172 | return IPMData(A, b, pb.objsense, c, c0, l, u) |
|
173 | - | end |
111 | 111 | # Lower-bound residual |
|
112 | 112 | # rl_j = l_j - (x_j - xl_j) if l_j ∈ R |
|
113 | 113 | # = 0 if l_j = -∞ |
|
114 | - | @. res.rl = (dat.l + pt.xl) - pt.x |
|
115 | - | res.rl .*= dat.lflag |
|
114 | + | @. res.rl = ((dat.l + pt.xl) - pt.x) * dat.lflag |
|
116 | 115 | ||
117 | 116 | # Upper-bound residual |
|
118 | 117 | # ru_j = u_j - (x_j + xu_j) if u_j ∈ R |
|
119 | 118 | # = 0 if u_j = +∞ |
|
120 | - | @. res.ru = dat.u - (pt.x + pt.xu) |
|
121 | - | res.ru .*= dat.uflag |
|
119 | + | @. res.ru = (dat.u - (pt.x + pt.xu)) * dat.uflag |
|
122 | 120 | ||
123 | 121 | # Dual residual |
|
124 | 122 | # rd = c - (A'y + zl - zu) |
188 | 186 | # Check for infeasibility certificates |
|
189 | 187 | if max( |
|
190 | 188 | norm(dat.A * pt.x, Inf), |
|
191 | - | norm((pt.x - pt.xl) .* dat.lflag, Inf), |
|
192 | - | norm((pt.x + pt.xu) .* dat.uflag, Inf) |
|
189 | + | norm((pt.x .- pt.xl) .* dat.lflag, Inf), |
|
190 | + | norm((pt.x .+ pt.xu) .* dat.uflag, Inf) |
|
193 | 191 | ) * (norm(dat.c, Inf) / max(1, norm(dat.b, Inf))) < - ϵi * dot(dat.c, pt.x) |
|
194 | 192 | # Dual infeasible, i.e., primal unbounded |
|
195 | 193 | mpc.primal_status = Sln_InfeasibilityCertificate |
197 | 195 | return nothing |
|
198 | 196 | end |
|
199 | 197 | ||
200 | - | δ = dat.A' * pt.y + (pt.zl .* dat.lflag) - (pt.zu .* dat.uflag) |
|
198 | + | δ = dat.A' * pt.y .+ (pt.zl .* dat.lflag) .- (pt.zu .* dat.uflag) |
|
201 | 199 | if norm(δ, Inf) * max( |
|
202 | 200 | norm(dat.l .* dat.lflag, Inf), |
|
203 | 201 | norm(dat.u .* dat.uflag, Inf), |
367 | 365 | # I. Recover positive primal-dual coordinates |
|
368 | 366 | δx = one(T) + max( |
|
369 | 367 | zero(T), |
|
370 | - | (-3 // 2) * minimum((pt.x - dat.l) .* dat.lflag), |
|
371 | - | (-3 // 2) * minimum((dat.u - pt.x) .* dat.uflag) |
|
368 | + | (-3 // 2) * minimum((pt.x .- dat.l) .* dat.lflag), |
|
369 | + | (-3 // 2) * minimum((dat.u .- pt.x) .* dat.uflag) |
|
372 | 370 | ) |
|
373 | - | pt.xl .= ((pt.x - dat.l) .+ δx) .* dat.lflag |
|
374 | - | pt.xu .= ((dat.u - pt.x) .+ δx) .* dat.uflag |
|
371 | + | @. pt.xl = ((pt.x - dat.l) + δx) * dat.lflag |
|
372 | + | @. pt.xu = ((dat.u - pt.x) + δx) * dat.uflag |
|
375 | 373 | ||
376 | 374 | z = dat.c - dat.A' * pt.y |
|
377 | 375 | #= |
385 | 383 | no | no | 0 | 0 | |
|
386 | 384 | ----+-----+--------+---------+ |
|
387 | 385 | =# |
|
388 | - | pt.zl .= ( z ./ (dat.lflag + dat.uflag)) .* dat.lflag |
|
389 | - | pt.zu .= (-z ./ (dat.lflag + dat.uflag)) .* dat.uflag |
|
386 | + | @. pt.zl = ( z / (dat.lflag + dat.uflag)) * dat.lflag |
|
387 | + | @. pt.zu = (-z / (dat.lflag + dat.uflag)) * dat.uflag |
|
390 | 388 | ||
391 | 389 | δz = one(T) + max(zero(T), (-3 // 2) * minimum(pt.zl), (-3 // 2) * minimum(pt.zu)) |
|
392 | 390 | pt.zl[dat.lflag] .+= δz |
58 | 58 | # Compute hx, hy, hz from first augmented system solve |
|
59 | 59 | hx = tzeros(Tv, n) |
|
60 | 60 | hy = tzeros(Tv, m) |
|
61 | - | ξ_ = dat.c - ((pt.zl ./ pt.xl) .* dat.l) .* dat.lflag - ((pt.zu ./ pt.xu) .* dat.u) .* dat.uflag |
|
61 | + | ξ_ = @. (dat.c - ((pt.zl / pt.xl) * dat.l) * dat.lflag - ((pt.zu / pt.xu) * dat.u) * dat.uflag) |
|
62 | 62 | KKT.solve!(hx, hy, hsd.kkt, dat.b, ξ_) |
|
63 | 63 | ||
64 | 64 | # Recover h0 = ρg + κ / τ - c'hx + b'hy - u'hz |
67 | 67 | h0 = ( |
|
68 | 68 | dot(dat.l .* dat.lflag, (dat.l .* θl) .* dat.lflag) |
|
69 | 69 | + dot(dat.u .* dat.uflag, (dat.u .* θu) .* dat.uflag) |
|
70 | - | - dot(c + (θl .* dat.l) .* dat.lflag + (θu .* dat.u) .* dat.uflag, hx) |
|
70 | + | - dot((@. (c + (θl * dat.l) * dat.lflag + (θu * dat.u) * dat.uflag)), hx) |
|
71 | 71 | + dot(b, hy) |
|
72 | 72 | + pt.κ / pt.τ |
|
73 | 73 | + hsd.regG |
77 | 77 | @timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0, |
|
78 | 78 | # Right-hand side of Newton system |
|
79 | 79 | res.rp, res.rl, res.ru, res.rd, res.rg, |
|
80 | - | -(pt.xl .* pt.zl) .* dat.lflag, |
|
81 | - | -(pt.xu .* pt.zu) .* dat.uflag, |
|
82 | - | -pt.τ * pt.κ |
|
80 | + | .-(pt.xl .* pt.zl) .* dat.lflag, |
|
81 | + | .-(pt.xu .* pt.zu) .* dat.uflag, |
|
82 | + | .-pt.τ * pt.κ |
|
83 | 83 | ) |
|
84 | 84 | ||
85 | 85 | # Step length for affine-scaling direction |
91 | 91 | @timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0, |
|
92 | 92 | # Right-hand side of Newton system |
|
93 | 93 | η .* res.rp, η .* res.rl, η .* res.ru, η .* res.rd, η * res.rg, |
|
94 | - | (-pt.xl .* pt.zl .+ γ * pt.μ .- Δ.xl .* Δ.zl) .* dat.lflag, |
|
95 | - | (-pt.xu .* pt.zu .+ γ * pt.μ .- Δ.xu .* Δ.zu) .* dat.uflag, |
|
94 | + | (.-pt.xl .* pt.zl .+ γ * pt.μ .- Δ.xl .* Δ.zl) .* dat.lflag, |
|
95 | + | (.-pt.xu .* pt.zu .+ γ * pt.μ .- Δ.xu .* Δ.zu) .* dat.uflag, |
|
96 | 96 | -pt.τ * pt.κ + γ * pt.μ - Δ.τ * Δ.κ |
|
97 | 97 | ) |
|
98 | 98 | α = max_step_length(pt, Δ) |
222 | 222 | ||
223 | 223 | @timeit hsd.timer "Δτ" Δ.τ = ( |
|
224 | 224 | ξg_ |
|
225 | - | + dot(dat.c |
|
226 | - | + ((pt.zl ./ pt.xl) .* dat.l) .* dat.lflag |
|
227 | - | + ((pt.zu ./ pt.xu) .* dat.u) .* dat.uflag |
|
225 | + | + dot((@. (dat.c |
|
226 | + | + ((pt.zl / pt.xl) * dat.l) * dat.lflag |
|
227 | + | + ((pt.zu / pt.xu) * dat.u) * dat.uflag)) |
|
228 | 228 | , Δ.x) |
|
229 | 229 | - dot(dat.b, Δ.y) |
|
230 | 230 | ) / h0 |
236 | 236 | ||
237 | 237 | # III. Recover Δxl, Δxu |
|
238 | 238 | @timeit hsd.timer "Δxl" begin |
|
239 | - | @. Δ.xl = -ξl + Δ.x - Δ.τ .* (dat.l .* dat.lflag) |
|
240 | - | Δ.xl .*= dat.lflag |
|
239 | + | @. Δ.xl = (-ξl + Δ.x - Δ.τ .* (dat.l .* dat.lflag)) * dat.lflag |
|
241 | 240 | end |
|
242 | 241 | @timeit hsd.timer "Δxu" begin |
|
243 | - | @. Δ.xu = ξu - Δ.x + Δ.τ .* (dat.u .* dat.uflag) |
|
244 | - | Δ.xu .*= dat.uflag |
|
242 | + | @. Δ.xu = ( ξu - Δ.x + Δ.τ .* (dat.u .* dat.uflag)) * dat.uflag |
|
245 | 243 | end |
|
246 | 244 | ||
247 | 245 | # IV. Recover Δzl, Δzu |
179 | 179 | ||
180 | 180 | # II. Recover Δxl, Δxu |
|
181 | 181 | @timeit mpc.timer "Δxl" begin |
|
182 | - | @. Δ.xl = -ξl + Δ.x |
|
183 | - | Δ.xl .*= dat.lflag |
|
182 | + | @. Δ.xl = (-ξl + Δ.x) * dat.lflag |
|
184 | 183 | end |
|
185 | 184 | @timeit mpc.timer "Δxu" begin |
|
186 | - | @. Δ.xu = ξu - Δ.x |
|
187 | - | Δ.xu .*= dat.uflag |
|
185 | + | @. Δ.xu = ( ξu - Δ.x) * dat.uflag |
|
188 | 186 | end |
|
189 | 187 | ||
190 | 188 | # III. Recover Δzl, Δzu |
259 | 257 | # Step length for affine-scaling direction |
|
260 | 258 | αp_aff, αd_aff = mpc.αp, mpc.αd |
|
261 | 259 | μₐ = ( |
|
262 | - | dot((pt.xl + αp_aff .* Δ.xl) .* dat.lflag, pt.zl + αd_aff .* Δ.zl) |
|
263 | - | + dot((pt.xu + αp_aff .* Δ.xu) .* dat.uflag, pt.zu + αd_aff .* Δ.zu) |
|
260 | + | dot((@. ((pt.xl + αp_aff * Δ.xl) * dat.lflag)), pt.zl .+ αd_aff .* Δ.zl) |
|
261 | + | + dot((@. ((pt.xu + αp_aff * Δ.xu) * dat.uflag)), pt.zu .+ αd_aff .* Δ.zu) |
|
264 | 262 | ) / pt.p |
|
265 | 263 | σ = clamp((μₐ / pt.μ)^3, sqrt(eps(T)), one(T) - sqrt(eps(T))) |
|
266 | 264 |
296 | 294 | αd_ = min(αd + δ, one(T)) |
|
297 | 295 | ||
298 | 296 | g = dot(pt.xl, pt.zl) + dot(pt.xu, pt.zu) |
|
299 | - | gₐ = dot((pt.xl + mpc.αp .* Δ.xl) .* dat.lflag, pt.zl + mpc.αd .* Δ.zl) + |
|
300 | - | dot((pt.xu + mpc.αp .* Δ.xu) .* dat.uflag, pt.zu + mpc.αd .* Δ.zu) |
|
297 | + | gₐ = dot((@. ((pt.xl + mpc.αp * Δ.xl) * dat.lflag)), pt.zl .+ mpc.αd .* Δ.zl) + |
|
298 | + | dot((@. ((pt.xu + mpc.αp * Δ.xu) * dat.uflag)), pt.zu .+ mpc.αd .* Δ.zu) |
|
301 | 299 | μ = (gₐ / g) * (gₐ / g) * (gₐ / pt.p) |
|
302 | 300 | ||
303 | 301 | # Newton RHS |
88 | 88 | # Lower-bound residual |
|
89 | 89 | # rl_j = τ*l_j - (x_j - xl_j) if l_j ∈ R |
|
90 | 90 | # = 0 if l_j = -∞ |
|
91 | - | @. res.rl = - pt.x + pt.xl + pt.τ * dat.l |
|
92 | - | res.rl .*= dat.lflag |
|
91 | + | @. res.rl = (- pt.x + pt.xl + pt.τ * dat.l) * dat.lflag |
|
93 | 92 | ||
94 | 93 | # Upper-bound residual |
|
95 | 94 | # ru_j = τ*u_j - (x_j + xu_j) if u_j ∈ R |
|
96 | 95 | # = 0 if u_j = +∞ |
|
97 | - | @. res.ru = - pt.x - pt.xu + pt.τ * dat.u |
|
98 | - | res.ru .*= dat.uflag |
|
96 | + | @. res.ru = (- pt.x - pt.xu + pt.τ * dat.u) * dat.uflag |
|
99 | 97 | ||
100 | 98 | # Dual residual |
|
101 | 99 | # rd = t*c - A'y - zl + zu |
173 | 171 | # Check for infeasibility certificates |
|
174 | 172 | if max( |
|
175 | 173 | norm(dat.A * pt.x, Inf), |
|
176 | - | norm((pt.x - pt.xl) .* dat.lflag, Inf), |
|
177 | - | norm((pt.x + pt.xu) .* dat.uflag, Inf) |
|
174 | + | norm((pt.x .- pt.xl) .* dat.lflag, Inf), |
|
175 | + | norm((pt.x .+ pt.xu) .* dat.uflag, Inf) |
|
178 | 176 | ) * (norm(dat.c, Inf) / max(1, norm(dat.b, Inf))) < - ϵi * dot(dat.c, pt.x) |
|
179 | 177 | # Dual infeasible, i.e., primal unbounded |
|
180 | 178 | hsd.primal_status = Sln_InfeasibilityCertificate |
182 | 180 | return nothing |
|
183 | 181 | end |
|
184 | 182 | ||
185 | - | δ = dat.A' * pt.y + (pt.zl .* dat.lflag) - (pt.zu .* dat.uflag) |
|
183 | + | δ = dat.A' * pt.y .+ (pt.zl .* dat.lflag) .- (pt.zu .* dat.uflag) |
|
186 | 184 | if norm(δ, Inf) * max( |
|
187 | 185 | norm(dat.l .* dat.lflag, Inf), |
|
188 | 186 | norm(dat.u .* dat.uflag, Inf), |
Files | Coverage |
---|---|
src | 89.02% |
Project Totals (43 files) | 89.02% |